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Abstract. 

This review paper presents an overview of the theoretical and experimental progress 
on the study of matter- wave dark solitons in atomic Bose-Einstein condensates. Upon 
introducing the general framework, we discuss the statics and dynamics of single and 
multiple matter-wave dark solitons in the quasi one-dimensional setting, in higher- 
dimensional settings, as well as in the dimensionality crossover regime. Special 
attention is paid to the connection between theoretical results, obtained by various 
analytical approaches, and relevant experimental observations. 
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1. Introduction 

A dark soliton is an envelope soliton having the form of a density dip with a phase 
jump across its density minimum. This localized nonlinear wave exists on top of a 
stable continuous wave (or extended finite-width) background. Dark solitons are the 
most fundamental nonlinear excitations of a universal model, the nonlinear Schrodinger 
(NLS) equation with a defocusing nonlinearity, and, as such, they have been studied in 
diverse branches of physics. Importantly, apart from a vast amount of literature devoted 
to relevant theoretical works, there exist many experimental results on dark solitons, 
including the observation of optical dark solitons, either as temporal pulses in optical 
fibers [UI2], or as spatial structures in bulk media and waveguides [SIS], the excitation 
of a non-propagating kink in a parametrically-driven shallow liquid [5], dark soliton 
standing waves in a discrete mechanical system [6], high-frequency dark solitons in thin 
magnetic films [TJ, dissipative dark solitons in a complex plasma [8J, and so on. 

Theoretical studies on dark solitons started as early as in 1971 [9] in the context 
of Bose-Einstein condensates (BECs). In particular, in Ref. [9], exact soliton solutions 
of the Gross-Pitaevskii (GP) equation (which is a variant of the NLS model) [10] were 
found and connected, in the small- amplitude limit, with the solitons of the Korteweg- 
de Vries (KdV) equation. Later, and shortly after the integration of the focusing 
NLS equation [TT], the defocusing NLS equation was also shown [12] to be completely 
integrable by means of the Inverse Scattering Transform (1ST) [13] ; this way, single as 
well as multiple dark soliton solutions of arbitrary amplitudes were found analytically. 
The 1ST approach allowed for an understanding of the formation of dark solitons [T4Tfl9"] . 
the interaction and collision between dark solitons [T2|[20] (see also Refs. [2lTf2"5] and [26] 
for relevant theoretical and experimental studies, respectively), and paved the way 
for the development of perturbation methods for investigating their dynamics in the 
presence of perturbations [25 1I271 - I32] . From a physical standpoint, dark solitons were 
mainly studied in the field of nonlinear optics — from which the term "dark" was coined. 
The first theoretical work in this context, namely the prediction of dark solitons in 
nonlinear optical fibers at the normal dispersion regime [33] , was subsequently followed 
by extensive studies of optical dark solitons [341135] . 

A new era for dark solitons started shortly after the realization of atomic BECs 
[361438] ; this achievement was awarded the Nobel prize in physics of 2001 [39|l4*0]. and 
has been recognized as one of the most fundamental recent developments in quantum 
and atomic physics over the last decades (see, e.g., the books jJTJHJ] f° r reviews). In an 
effort to understand the properties of this exciting state of matter, there has been much 
interest in the macroscopic nonlinear excitations of BECs (see reviews in Refs. [4*3ll4*4"] ) . 
In that regard, the so-called matter-wave dark solitons, were among the first purely 
nonlinear states that were experimentally observed in BECs [4"5H4*9] . 

The interest on matter- wave dark solitons is not surprising due to a series of reasons. 
First of all, for harmonically confined BECs, these structures are the nonlinear analogues 
of the excited states of a "prototype" quantum system [SUET], namely the quantum 
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harmonic oscillator |52j. On the other hand, the topological nature of matter- wave dark 
solitons (due to the phase jump at their density minimum) renders them a "degenerate" , 
one-dimensional (ID) analogue of vortices, which are of paramount importance in diverse 
branches of physics [53]. Additionally, and perhaps more importantly, matter- wave 
dark solitons are — similarly to vortices [53H56] — quite fundamental structures arising 
spontaneously upon crossing the BEC phase-transition [57| l58]. with properties which 
may be used as diagnostic tools probing the rich physics of a purely quantum system 
(BEC) at the mesoscale [59]. Finally, as concerns applications, it has been proposed 
that the dark soliton position can be used to monitor the phase acquired in an atomic 
matter- wave interferometer in the nonlinear regime [601I6T] (see also relevant experiments 
of Refs. [62,63j devoted to atom-chip interferometry of BECs). 

The early matter-wave dark soliton experiments, as well as previous works on dark 
solitons in optics, inspired many theoretical efforts towards a better understanding of 
the stability, as well as the static and dynamical properties of matter-wave dark solitons. 
Thus, it is probably not surprising that a new series of experimental results from various 
groups have appeared [MHTT] . while still other experiments - - not directly related 
to dark solitons — reported observation of these structures [62],[63j[72]. These new, 
very recent, experimental results were obtained with an unprecedented control over the 
condensate and the solitons as compared to the earlier soliton experiments. Thus, these 
"new age" experiments were able not only to experimentally verify various theoretical 
predictions, but also to open new exciting possibilities. Given this emerging interest, 
and how new experiments in BEC physics inspire novel ideas — both in theory and in 
experiments — new exciting results are expected to appear. 

The present paper aims to provide an overview of the theoretical and experimental 
progress on the study of dark solitons in atomic BECs. The fact that there are 
many similarities between optical and matter- wave dark solitons [73], while there exist 
excellent reviews on both types of dark solitons (see Ref. [31] for optical dark solitons 
and Ch. 4 in Ref. [13] for matter- wave dark solitons), provides some restrictions in 
the article: first, the space limitations of the article, will not allow for an all-inclusive 
presentation; in that regard, important entities — relevant to dark solitons — such as 
vortices [53l[74j[75] an d vortex rings [SEE?] will only be discussed briefly. In fact, this 
review (which obviously entails a "personalized" perspective on dark solitons) will cover 
the basic theory emphasizing, in particular, to the connection between the theoretical 
results and experimental observations; this way, in most cases, theoretical discussion 
will be immediately followed by a presentation of pertinent experimental results. In 
that regard, it is also relevant to note that our theoretical approach will basically be 
based on the mean-field theory: as will be shown, the latter can be used as a basis of 
understanding of most effects and experimental findings related to matter-wave dark 
solitons; this way, thermal and quantum effects — which may be particularly relevant 
and important in certain cases — will only be briefly covered. Following the above 
limitations, the structure of the manuscript will be as follows. 

Section [2] is devoted to the mean-field description of BECs. Particularly, we first 
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present the GP equation and discuss its connection with the respective full quantum 
many-body problem. Next, we present the ground state of the condensate and discuss 
how its small-amplitude excitations can be studied by means of the Bogoliubov-de 
Gennes (BdG) equations. Lower-dimensional versions of the GP model, pertinent to 
highly anisotropic trapping potentials, are also discussed; this way, depending on the 
shape of the trap, we start from purely three-dimensional (3D) BECs and introduce 
elongated (alias "cigar-shaped") BECs, quasi one-dimensional (ID) BECs and quasi 
two-dimensional (2D) (alias "disk-shaped") ones, as well as discuss cases relevant to 
the dimensionality crossover regimes. The topics of strongly-interacting Bose gases, and 
their relevant mean-field description, are also briefly covered. 

Section [3] provides the theoretical basis for the study of matter- wave dark solitons. 
Specifically, first we present the completely integrable ID NLS equation, its basic 
properties and the dark soliton solutions. Relevant mathematical tools, such as Inverse 
Scattering Transform (1ST), the renormalization of the integrals of motion of dark 
solitons and the small-amplitude approximation — leading to the connection of matter- 
wave dark solitons to Korteweg-de Vries (KdV) solitons — are discussed. Furthermore, 
the generation of matter-wave dark solitons by means of the phase-, density- and 
quantum-state-engineering methods are also presented. We also provide the multiple- 
dark soliton solutions of the NLS equation, and discuss their interactions and collisions. 

Section H] deals with matter-wave dark solitons in quasi-lD Bose gases. Particularly, 
we first discuss the adiabatic dynamics of dark solitons in the presence of the harmonic 
trap by means of different analytical techniques; these include the Hamiltonian and 
Lagrangian approaches of the perturbation theory, the Landau dynamics and the small- 
amplitude approximation approaches. Next, a connection between the stability, statics 
and dynamics of dark solitons is presented, relying on a study of the Bogoliubov 
spectrum of single- and multiple-dark solitons and the role of the pertinent anomalous 
modes. Non-adiabatic effects, namely emission of radiation of solitons in the form 
of sound waves, as well as rigorous results concerning the persistence and stability of 
matter-wave dark solitons, are also discussed. 

Section [5] studies matter-wave dark solitons in higher-dimensional settings. 
Considering, at first, the case of purely 2D or 3D geometries, the transverse (alias 
"snaking") instability of rectilinear dark solitons, and the concomitant soliton decay 
into vortex pairs or vortex rings, is presented. The theme of matter-wave dark solitons 
of radial symmetry, namely ring dark solitons and spherical shell solitons, is also covered. 
Furthermore, we present results concerning the stability of dark solitons in cigar-shaped 
(3D) BECs, and in BECs in the dimensionality crossover regime from 3D to ID; in the 
latter experimentally relevant setting, both single- and multiple- dark soliton statics and 
dynamics are analyzed. 

In Section El we discuss various experimentally relevant settings and parameter 
regimes for matter-wave dark solitons. In particular, we first present results concerning 
matter- wave dark solitons in multi-component (pseudo-spinor and spinor) BECs. Next, 
we discuss how matter-wave interference and the breakdown of BEC superfluidity are 



6 



connected to the generation of matter-wave dark solitons. We continue by referring 
to matter-wave dark solitons in periodic potentials, namely optical lattices (OLs) and 
superlattices, and conclude this Section by discussing the statics and dynamics of dark 
solitons at finite temperatures. 

Finally, in Section [7] we briefly summarize our conclusions and discuss future 
challenges. 

2. Mean-field description of Bose-Einstein condensates 

Bose-Einstein condensation of dilute atomic gases is an unambiguous manifestation 
of a macroscopic quantum state in a many-body system. As such, this phenomenon 
has triggered an enormous amount of experimental and theoretical work [41, 42J. 
Importantly, this field is intimately connected with branches of physics such as 
superfluidity, superconductivity, lasers, coherent optics, nonlinear optics, and physics 
of nonlinear waves. Many of the common elements between BEC physics and the above 
areas, and in particular optics, rely on the existence of macroscopic coherence in the 
many-body state of the system. From a theoretical standpoint, this can be understood 
by the fact that many effects related to BEC physics can be described by a mean- 
field model, namely the Gross-Pitaevskii (GP) equation [10]. The latter is a partial 
differential equation (PDE) of the NLS type, which plays a key role — among other 
fields — in nonlinear optics [35]. Thus, BEC physics is closely connected to nonlinear 
optics (and the physics of nonlinear waves), with vortices and solitons being perhaps the 
most prominent examples of common nonlinear structures arising in these areas [431144] . 

Below we will briefly discuss the theoretical background for the description of BECs. 
We emphasize, in particular, lowest-order mean-field theory, as this can be used as a 
basis to understand the nonlinear dynamics of matter-wave dark solitons. Interesting 
effects naturally arise beyond the GP mean-field, both due to thermal and quantum 
fluctuations. Such effects become particularly relevant in extremely tightly confining 
geometries, or when the Bose-Einstein condensation transition region is approached. 

2.1. The Gross-Pitaevskii equation. 

In order to describe theoretically the statics and dynamics of BECs a quantum many- 
body approach is required [411142] (see also Ref. [78] for a recent review on the many-body 
aspects of BECs). Particularly, a sufficiently dilute ultracold atomic gas, composed by N 
interacting bosons of mass m confined by an external potential V^ xt (r), can be described 
by the many-body Hamiltonian; the latter can be expressed, in second quantization 
form, through the boson annihilation and creation field operators, \T/(r, t) and ^'(r, t) 
(which create and annihilate a particle at the position r) namely, 
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(1) 
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where Hq = — (h 2 /2m)V 2 + Kxt( r ) is the single-particle operator and V(r — r') is the 
two-body interatomic potential. Apparently, the underlying full many-body problem is 
very difficult to be treated analytically (or even numerically) as N increases and, thus, 
for convenience, a mean-field approach can be adopted. The mean-field approach is 
based on the separation of the condensate contribution from the boson field operator as 
follows 



(r, *) = (tt(r, t)) + #'(r, t) = *(r, *) + #'(r, t). (2) 

In the above expression, the expectation value of the field operator, (ijf(r,t)) = ^(r, t), 
is known as the macroscopic wave function of the condensate, while \&'(r, t) describes the 
non-condensate part, which accounts for quantum and thermal fluctuations. Considering 
the case of a dilute ultracold gas with binary collisions at low energy, characterized by 
the s-wave scattering length a, the interatomic potential can be replaced by an effective 
delta-function interaction potential, V(r' — r) = gS(r' — r) [U1H2~] . with the coupling 
constant g given by g = 4irh 2 a/m. Under these assumptions, a nontrivial zeroth-order 
theory for the BEC wave function can be obtained by means of the Heisenberg evolution 
equation ih(d4f /dt) = [V, H], upon replacing the field operator ^ with the classical field 
i.e., ignoring the quantum and thermal fluctuations described by if>'(r',t). Such a 
consideration leads to the Gross-Pitaevskii (GP) equation [10] . which has the form: 



ihd t ^{r,t) 



V 2 

2m 



Kxt(r) + <?I*M)P 



*(r,t). 



In the above equation, ^(r, t) is normalized to the number of atoms N, namely, 



N 



\^(r,t)\ 2 dr, 



(3) 



(4) 



and the nonlinearity (which is obviously introduced by interatomic interactions) is 
characterized by the s-wave scattering length a, which is a > or a < for repulsive 
or attractive interatomic interactions, respectively. Notice that Eq. can be written 
in canonical form, ihdt^ = 5E/5^* (with * denoting complex conjugate), where the 
dynamically conserved energy functional E is given by 



E 



dr 



^iwi 2 +Kxti^r+^i 4 

2m 2 



(5) 



with the three terms in the right-hand side representing, respectively, the kinetic energy, 
the potential energy and the interaction energy. 

A time-independent version of the GP equation can be obtained upon expressing the 
BEC wave function as \I/(r, t) = ^o(r) exp(—i/j,t/H), where \i = dE/dN is the chemical 
potential. This way, Eq. ([3]) yields the following equation for the stationary state \l/o : 

h 2 



2m 



V 2 + V ext (r)+g\y \ 2 (r) 



tf (r) = ^tfo(r) 



(6) 
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2.2. The mean-field approach vs. the many-body quantum mechanical problem. 

Although the GP equation is known from the early 60's jlQ] . it was only recently shown 
that it can be derived rigorously from a self-consistent treatment of the respective many- 
body quantum mechanical problem [80] . In particular, in Ref . [80] — which dealt with 
the stationary GP Eq. ([6]) — it was proved that the GP energy functional describes 
correctly the energy and the particle density of a trapped Bose gas to the leading-order 
in the small parameter n\a\ 3 f , where n is the average density of the gas. The above 
results were proved in the limit where the number of particles iV — > oo and the scattering 
length a — > 0, such that the product Na stays constant. Importantly, although Ref. [80] 
referred to the full three-dimensional (3D) Bose gas, extensions of this work for lower- 
dimensional settings were also reported (see the review [81] and references therein). 

The starting point of the analysis of Ref. [80] is the effective Hamiltonian of iV 
identical bosons, which can be expressed (in units so that h~ = 2m = 1) as follows: 



where v (|r|) is a general interaction potential assumed to be spherically symmetric and 
decaying faster than |r|~ 3 at infinity. Then, assuming that the quantum- mechanical 
ground-state energy of the Hamiltonian (GO) is Eqm(N, a) (here N is the number of 
particles and a is the dimensionless two-body scattering length), the main theorem 
proved in Ref. [80J is the following. The GP energy is the dilute limit of the quantum- 
mechanical energy: 



where Eqj>(N, a) is the energy of a solution of the stationary GP Eq. (jH]) (in units such 
that h = 2m = 1), and the convergence is uniform on bounded intervals of 5i. 

The above results (as well as the ones in Ref. [81] ) were proved for stationary 
solutions of the GP equation, and, in particular, for the ground state solution. More 
recently, the time-dependent GP Eq. ([3]) was also analyzed within a similar asymptotic 
analysis in Ref. [82]. In this work, it was proved that the limit points of the fc-particle 
density matrices of ^N,t (which is the solution of the iV-particle Schrodinger equation) 
satisfy asymptotically the GP equation (and the associated hierarchy of equations) with 
a coupling constant given by J v(x)dx, where v(x) describes the interaction potential. 

These rigorous results, as well as a large number of experimental results related to 
the physics of BECs, indicate that (under certain conditions) the GP equation is a good 
starting point for analyzing the statics and dynamics of BECs. 

f The condition n\a\ 3 -C 1, which is also required for the derivation of the GP Eq. ([3]), implies that 
the Bose-gas is "dilute" or "weakly-interacting"; typically, in BEC experiments, n\a\ 3 < 10 -3 [42] . 



N 




(7) 





9 



2.3. Ground state and excitations of the condensate. 



Let us now consider a condensate confined in a harmonic external potential, namely, 



Kxt(r) = -m(ulx 2 + uj 2 y y 2 + u 2 z z 2 ), 



(9) 



where u 7 



UJ 



and uj z are the (generally different) trap frequencies along the three 



directions. In this setting, and in the case of repulsive interatomic interactions 
(a > 0) and sufficiently large number of atoms N, Eq. (|6]) can be used to determine 
analytically the ground state of the system. In particular, in the asymptotic limit of 
Na/a^o ^> 1 (where aho = yhf (rnuho) is the harmonic oscillator length associated 
with the geometrical average Uho = (ux^yUz) of the trap frequencies), it is expected 
that the atoms are pushed towards the rims of the condensate, resulting in slow spatial 
variations of the density profile n(r) = |\l/o( r )| 2 - Thus, the latter can be obtained as an 
algebraic solution stemming from Eq. when neglecting the kinetic energy term - 
the so-called Thomas-Fermi (TF) limit pTIT - H3] : 

n(r) = g- 1 [ f i-V cxt (r)}, (10) 

in the region where \i > V ext (r), and n = outside, and the value of \i being determined 
by the normalization condition [cf. Eq. (j4])]. Notice that the TF approximation becomes 
increasingly accurate for large values of \i. 

Small- amplitude excitations of the BEC can be studied upon linearizing Eq. 
around the ground state. Particularly, we consider small perturbations of this state, i.e., 



tt(r,t) 



3 — i/tt//i 



*o(r) + £ ( 



//fir)e-^* + u;(r)e^*) 



:ir: 



where Uj, Vj are the components of the linear response of the BEC to the external 
perturbations that oscillate at frequencies ±uij [the latter are (generally complex) 
eigenfrequencies] . Substituting Eq. (TTTT) into Eq. and keeping only the linear terms 
in Uj and Vj, we obtain the so-called Bogoliubov-de Gennes (BdG) equations: 

H -fi + 2g |^ | 2 (r)l u,(r) + g *g(r)uj(r) 



hujjUjir) 



H -ti + 2g\y \ 2 {r) 



g% 2 [ 



T)Uj[T) 



-flUJjVjiv) 



(12) 



where Hq = — (h 2 /2m)'V 2 + V ext (r) is the single-particle Hamiltonian. Importantly, 
these equations can also be used, apart from the ground state, for any other stationary 
state (including, e.g., solitons) with the function being modified accordingly. In such 
a general context, the BdG equations provide the eigenfrequencies u = uj r + iui and 
the amplitudes Uj and Vj of the normal modes of the system. Note that due to the 
Hamiltonian nature of the system, if u is an eigenfrequency of the Bogoliubov spectrum, 
so are —u, u* and — uj* . In the case of stable configurations with ui = 0, the solution 
of BdG equations with frequency u represent the same physical oscillation with the 
solution with frequency —u [12]. 



10 



In the case of a homogeneous gas (V cxt (r) = 0) characterized by a constant density 
n = l^ol 2 ; the amplitudes Uj and Vj in the BdG equations are plane waves, ~ exp(zk-r), 
of wave vector k. Then, Eqs. ( TT2l) lead to the dispersion relation, 

w 2 =(^)(^r + H- (13) 

In the case of repulsive interatomic interactions (g > 0), Eq. ( f!3|) indicates that small- 
amplitude harmonic excitations of the stationary state 

$ = y/n^exp(-ifit/h), (14) 

with fi — no, are always stable since cjj = for every k. Thus, this state is not subject 
to the modulational instability (see, e.g., Ref. [S3] and references therein). This fact is 
important, as the wave function of Eq. (I14p can serve as a stable background (alias 
"pedestal"), on top of which strongly nonlinear localized excitations may be formed; 
such excitations may be, e.g., matter- wave dark solitons which are of particular interest 
in this work. Notice that the above mentioned small-amplitude harmonic excitations 
are in fact sound waves, characterized by the phonon dispersion relation u = |k|c s [see 
Eq. (fl3|) for small momenta Kk], where 

c s = ^gn Q /m, (15) 

is the speed of sound. We should note in passing that in the case of attractive interatomic 
interactions (g < 0), the speed of sound becomes imaginary, which indicates that long 
wavelength perturbations grow or decay exponentially in time. Thus, the stationary 
state of Eq. (fl4"]) is subject to the modulational instability, which is responsible for 
the formation of matter-wave bright solitons [841486] in attractive BECs (see also the 
reviews [43]I111[E3],[EZ] and references therein). 

2.4- Lower- dimensional condensates and relevant mean-field models. 

Let us consider again a condensate confined in the harmonic trap of Eq. (Q. In this 
case, the trap frequencies set characteristic length scales for the spatial size of the 
condensate through the harmonic oscillator lengths aj = (h/muj) 1 ^ 2 (j G {x, y, z}). 
Another important length scale, introduced by the effective mean-field nonlinearity, is 
the so-called healing length defined as £ = (87m a) _1 ^ 2 (with n being the maximum 
condensate density). The healing length, being the scale over which the BEC wave 
function "heals" over defects, sets the spatial widths of nonlinear excitations, such as 
matter-wave dark solitons. 

Based on the above, as well as the form of the ground state [cf. Eq. f fTOj) ]. it is clear 
that the shape of the BEC is controlled by the relative values of the trap frequencies. 
For example, if u x = u y = u>± ~ u z (i.e., for an isotropic trap), the BEC is almost 
spherical, while for u z < u± (i.e., for an anisotropic trap) the BEC is "cigar shaped". 
It is clear that such a cigar-shaped BEC (a) may be a purely 3D object, (b) acquire an 
almost ID character (for strongly anisotropic traps with u z u± and /i hw±), or (c) 
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being in the so-called dimensionality crossover regime from 3D to ID. These regimes 
can be described by the dimensionless parameter [88J, 

d = NQ— , (16) 
a± 

where Q = u z /u± is the so-called "aspect ratio" of the trap. Particularly, if the 
dimensionality parameter is d ^> 1, the BEC locally retains its original 3D character 
(although it may have an elongated, quasi- ID shape) and its ground state can be 
described by the TF approximation in all directions. On the other hand, if d <C 1, 
excited states along the transverse direction are not energetically accessible and the 
BEC is effectively ID. Apparently, this regime is extremely useful for an analytical 
study of matter-wave dark solitons. Finally, if d ~ 1, the BEC is in the crossover regime 
between ID and 3D, which is particularly relevant as recent matter-wave dark soliton 
experiments have been conducted in this regime [69|I7T]. 

Let us now discuss in more detail lower-dimensional mean-field models describing 
cigar-shaped BECs. First, we consider the quasi-lD regime (d <C 1) characterized by 
an extremely tight transverse confinement. In this case, following Refs. [50,89,90j, the 
BEC wave function is separated into transverse and longitudinal components, namely 
\I/(r, t) = §(r;t)ip(z,t). Then, the transverse component <E>(r;i) is described by the 
Gaussian ground state of the transverse harmonic oscillator (and, thus, the transverse 
width of the condensate is set by the transverse harmonic oscillator length a±), while 
the longitudinal wave function ip(z,t) obeys the following effectively ID GP equation: 

-^d? + V(z)+g 1D mz,t)\ 



ihd t ip(z } t) 



(17) 



2m 

where the effective ID coupling constant is given by gm = g/2na\ = 2ahw± and 
V(z) = (1 / '2)mu 2 z 2 . Notice that in the case under consideration, if the additional 
condition [(N/y/U)(a/a±)] 1 / 3 ^> 1 is fulfilled, then the longitudinal condensate density 
n(z, t) = \tp(z, t)\ 2 can be described by the TF approximation — see Eq. f lTUj) with \x 
now being the ID chemical potential (and g — >■ g\o) [88J. Following the terminology of 
Ref. [69], this regime will hereafter be referred to as the TF-1D regime. 

Next, let us consider the effect of the deviation from one-dimensionality on the 
longitudinal condensate dynamics. In this case, the wave function can be factorized 
as before, but with the transverse component $ assumed to depend also on the 
longitudinal variable z (and time t) [9~TH9~1] . Physically speaking, it is expected that 
the transverse direction will no longer be occupied by the ground state, but $ would 
still be approximated by a Gaussian function with a width w = w(z,t) that can be 
treated as a variational parameter [92TI91] . This way, it is possible to employ different 
variational approaches and derive the following NLS equation for the longitudinal wave 
function, 



h 2 d 2 



+ V{z) + f(n) 



2m dz 2 

The nonlinearity function f(n) in Eq. ([IB"]) depends on the longitudinal density n(z,t) 
and may take different forms. Particularly, in Ref. [92] (where variational equations 
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related to the minimization of the action functional were used) f(n) is found to be: 

f(n) = - \ . = + VI + 2an , 19 

and the respective NLS equation is known as the non-polynomial Schrodinger equation 
(NPSE). On the other hand, in Refs. [931194"] (where variational equations related to the 
minimization of the transverse chemical potential were used) the result for f(n) is: 

f(n) = hou ± Vl + Aan. (20) 

Since, as explained above, the derivation of the mean-field models with the nonlinearity 
functions in Eqs. (fT9|) and (|20|) is based on different approaches, these nonlinearity 
functions are quite different. Nevertheless, they can be "reconciled" in the weakly- 
interacting limit of an C 1: in this case, the width of the transverse wave function 
becomes w = a± and Eq. ( fl8l) — with either the nonlinearity function of Eq. ( fl9l) or 
that of Eq. (120]) — is reduced to the ID GP model of Eq. (TT7T) . 

The above effective ID models predict accurately ground state properties of quasi- 
1D condensates, such as the chemical potential, the axial density profile, the speed of 
sound, collective oscillations, and others. Importantly, these models are particularly 
useful in the dimensionality crossover regime, where they describe the axial dynamics 
of cigar-shaped BECs in a very good approximation to the 3D Gross-Pitaevskii (GP) 
equation (see, e.g., theoretical work related to matter- wave dark solitons in Ref. [95] 
and relevant experimental results in Ref. [69J). 

On the other hand, extremely weak deviations from one-dimensionality can also 
be treated by means of a rather simple non-cubic nonlinearity that can be obtained by 
Taylor expanding f(n), namely: 

f(n) = gin- g 2 n 2 , (21) 

where g\ = gm and g 2 depends on the form of f{n). In this case, Eq. (fl8l) becomes 
a cubic-quintic NLS (cqNLS) equation. This model was derived self-consistently 
in Ref. [9T] , where dynamics of matter- wave dark solitons in elongated BECs was 
considered; there, the coefficient g 2 was found to be equal to g 2 = 241n(4/3)a 2 ^wj_. 

Here, it is worth mentioning that the quintic term in the cqNLS equation may have 
a different physical interpretation, namely to describe three-body interactions, regardless 
of the dimensionality of the system. In this case, the coefficients gin and g 2 in Eq. fTST]) 
are generally complex, with the imaginary parts describing inelastic two- and three-body 
collisions, respectively |96j. As concerns the rate of the three-body collision process, it is 
given by (dn/dt) = —Ln 3 [41], where L is the loss rate (which is of order of 10~ 27 -10 -30 
cm 6 s _1 for various species of alkali atoms [97J). Accordingly, the decrease of the density 
is accounted for by the term — (L/2)\ip\ 4: ip in the time dependent GP equation, i.e., to 
the quintic term in the cqNLS equation. 

It is also relevant to note that the NLS Eq. ( TT8|) has also been used as a mean-field 
model describing strongly-interacting ID Bose gases and, particularly, the so-called 
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Tonks- Girardeau gas of impenetrable bosons [HE] (see also Refs. [221 100J for recent 
experimental observations). In this case, the function f(n) takes the form [101] . 

/(«) = ^" 2 , (22) 

and, thus, Eq. ( fl8i) becomes a quintic NLS equation. Although the applicability of 
this equation has been criticized (as in certain regimes it fails to predict correctly the 
coherence properties of the strongly- interacting ID Bose gases [102] ). the corresponding 
hydrodynamic equations for the density n and the phase (p arising from the quintic NLS 
equation under the Madelung transformation ip = y / nexp(z<^) are well-documented in 
the context of the local density approximation [103] . Additionally, it should be noticed 
that the time-independent version of the quintic NLS equation has been rigorously 
derived from the many-body Schrodinger equation [104] . 

We finally mention that another lower-dimensional version of the fully 3D GP 
equation can be derived for "disk-shaped" (alias "pancake") condensates confined in 
strongly anisotropic traps with u± <C uj z and fi <C Hcu z . In such a case, a procedure 
similar to the one used for the derivation of Eq. (ITT|) leads to the following (2 + 1)- 
dimensional NLS equation: 

" h 2 

-—Vl + V{r)+g 2D \^{x,y,t)[ 



ihd t ip(x,y,t) 



il)(x,y,t), (23) 



2m 

where r 2 = x 2 + y 2 , V 2 ± = d 2 + d 2 , the effectively 2D coupling constant is given by 
92D — g/VZftO'z — 2y / 27r aa z hco z , while the potential is given by V(r) = (l/2)mu}r 2 . It 
should also be noticed that other effective 2D mean-field models (involving systems of 
coupled 2D equations [92J or 2D GP equations with generalized nonlinearities [92] 194"]) 
have also been proposed for the study of the transverse dynamics of disk-shaped BECs. 

The above models will be used below to investigate the static and dynamical 
properties of matter-wave dark solitons arising in the respective settings. 

3. General background for the study of matter-wave dark solitons 

3.1. NLS equation and dark soliton solutions. 

We start by considering the case of a quasi-lD condensate described by Eq. (IPT1) . The 
latter, can be expressed in the following dimensionless form: 

1 



id t ip(z,t) 



-~d 2 + V(z) + mzM 



iP(z,t), (24) 



where the density \ip\ 2 , length, time and energy are respectively measured in units of 
2a, a_i_, wj 1 and hu±, while the potential V(z) is given by 

V(z) = ^V. (25) 

In the case under consideration, the normalized trap strength (aspect ratio) is Q 1 
and, thus, as a first step in our analysis, the potential V(z) is ignored, f In such a 

f Note that in the limit of z — > ±oo this approximation always breaks down. 
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case, the condensate is homogeneous and can be described by the completely integrable 
defocusing NLS equation [12] (see also the review [34J): 

1 



id t ip(z,t) 



-ld 2 z + \iP(z,t)(' 



1>(z,t). (26) 



This equation possesses an infinite number of conserved quantities (integrals of motion); 
the lowest-order ones are the number of particles N, the momentum P, and the energy 
E, respectively given by: 

/— oo 
\*P\ 2 dz, (27) 
-oo 
• /■- oo 

P = - / (*pd z r - dz, (28) 



2 



oo 
— oo 



E = \l (m\ 2 + \^)dz. (29) 

It is also noted that the NLS Eq. f j26|) can be obtained by the Euler-Lagrange equation 
5£/5ip* = d t (dd tt p*C) + d z (dd z ^*C)d^*C=0, where the Lagrangian density C is given by: 

c = \ - ro^) - \ o<^i 2 + m 4 ) . (30) 

The simplest nontrivial solution of Eq. ( 126]) is a plane wave of wave number k and 
frequency u>, namely, 

= \f^ey^>[i{k z ~ + ®o)], w = -k 2 — /i, (31) 

where the constant BEC density n sets the chemical potential, i.e., n = fi and 9 Q is 
an arbitrary constant phase. This solution, which is reduced to the stationary state 
of Eq. fTH|) for k = 0, is also modulationally stable as can be confirmed by a simple 
stability analysis (see, e.g., Refs. (3UIH3]). For small densities, no <C 1, the above plane 
wave satisfies the linear Schrodinger equation, idtip + ^d 2 ip = 0, and the pertinent linear 
wave solutions of the NLS equation are characterized by the dispersion relation oj = \k 2 . 
Notice that if the system is characterized by a length L, then the integrals of motion 
for the stationary solution in Eq. fj3T|) take the values: 

N = 2n L, P = kn L, E = ^(k 2 - n )n L. (32) 

The NLS equation admits nontrivial solutions, in the form of dark solitons, which 
can be regarded as strongly nonlinear excitations of the plane wave solution fl3T|) . In 
the most general case of a moving background [k ^ in Eq. ( 131]) ]. a single dark soliton 
solution may be expressed as [12J, 

ip(Zj t) = y/no (£?tanh£ + iA) exp[i(kz — cut + 9 )}, (33) 

where ( = ^JuqB [z — z (t)}; here, z (t) = vt + z Q is the soliton center, z Q is an arbitrary 
real constant representing the initial location of the dark soliton, v is the relative velocity 
between the soliton and the background given by v — A^/hq + k, the frequency uj is 
provided by the dispersion relation of the background plane wave, uj = (l/2)fc 2 + n 
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1 



-1 



A(j>=71 



-v=0 



-v=0.6c 



Figure 1. (color online) Examples of the density (top panel) and phase (bottom 
panel) of a black (blue line) and a gray (green line) soliton on top of a background 
with density no — 1. The black soliton's parameters are A = and B = 1, i.e., u = 0, 
"mm = and A0 = 7r. The gray soliton's parameters are A = 0.6 and B = 0.8, i.e., 
v = 0.6c s (here c s = = 1), n m in = n (l — B 2 ) = n A 2 = 0.36, and Acj) = 0.3i7r. 



[cf. Eq. ( 13~TT) ] f , and, finally, the parameters A and 5 are connected through the 
equation A 2 + B 2 = 1. In some cases it is convenient to use one parameter instead of 
two and, thus, one may introduce 

A = sin0, £> = cos0, (34) 

where <fi is the so-called "soliton phase angle" {\<p\ < tc/2). Notice that although the 
asymptotics of the dark soliton solution ( 1331) coincide with the ones of Eq. ( l3il . the 
plane waves at z — > ±oo have different phases; as a result, there exists a nontrivial 
phase jump A0 across the dark soliton, given by: 

B\ _ vr" 

A) 2 

Note that, hereafter, we will consider the simpler case where the background of matter- 
wave dark solitons is at rest, i.e., k = 0; then, the frequency u> actually plays the role 
of the normalized chemical potential, namely u = fi = n , which is determined by the 
number of atoms of the condensate. 

The soliton phase angle describes also the darkness of the soliton, namely, 

|^| 2 = n o (l-cos 2 0sech 2 C). (36) 

f Here, this dispersion relation implies that lo > k 2 and, thus, the allowable region in the (k, to) plane 
for dark solitons is located above the parabola u = ^k 2 corresponding to the linear waves. 



A0 



tan 



^tan- 1 !^)- (35) 
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This way, the cases = and < <fi < tt/2 correspond to the so-called black and 
gray solitons, respectively. The amplitude and velocity of the dark soliton are given (for 



is characterized by a zero velocity, v = (and, thus, it is also called stationary kink), 
while the gray soliton moves with a finite velocity v ^ 0. Examples of the forms of a 
black and a gray soliton are illustrated in Fig. [TJ 

In the limiting case of a very shallow (small- amplitude) dark soliton with cos <fi <C 1, 
the soliton velocity is close to the speed of sound which, in our units, is given by: 



The speed of sound is, therefore, the maximum possible velocity of a dark soliton which, 
generally, always travels with a velocity less than the speed of sound. We finally note 
that the dark soliton solution ( 1331) has two independent parameters (for k = 0), one for 
the background, no, and one for the soliton, 0, while there is also a freedom (translational 
invariance) in selecting the initial location of the dark soliton z f . 

In the case of a condensate confined in a harmonic trap [cf. Eq. (jHJ)], the background 
of the dark soliton is, in fact, of finite extent, being the ground state of the BEC [which 
may be approximated by the Thomas- Fermi cloud, cf. Eq. (flOf)]. For example, in the 
quasi- ID setting of the ID GP Eq. (124)) with the harmonic potential in Eq. (125]) . the 
"composite" wave function (describing both the background and the soliton) can be 
approximated as ip(z,t) = <3>(z) exp(— ifit)^As{z, t), where $(z) is the TF background 
and ip^s(z,t) is the dark soliton wave function of Eq. (l33l . which satisfies the ID GP 
equation for V(z) = 0. 

3.2. Dark solitons and the Inverse Scattering Transform. 

The single dark soliton solution of the NLS Eq. ( )26l) presented in the previous Section, 
as well as multiple dark soliton solutions (see Sec. 13.61 below), can be derived by means 
of the Inverse Scattering Transform (1ST) |12j . A basic step of this approach is the 
solution of the Zakharov-Shabat (ZS) eigenvalue problem, with eigenvalue A, for the 
auxiliary two-component eigenfunction U = (ui,U2) T , namely, 



with the boundary conditions i/;(z, 0) — > ^fr^o, for z — > +oo, and ip(z, 0) — > i/noex-p^iO), 
for z — > — oo. Here, y/no is the amplitude of the background wave function and 9 
is a constant phase. Since the operator L is self-adjoint, the ZS eigenvalue problem 
possesses real discrete eigenvalues Xj, with magnitudes |A 3 -| < ^JnQ. Importantly, each 
real discrete eigenvalue Xj = ^/n^sin^,,- corresponds to a dark soliton of depth ^/hq cos 4>j 

f Recall that the underlying model, namely the completely integrable NLS equation, has infinitely 
many symmetries, including translational and Galilean invariance. 



k = 0) by y^cos^ and y^sin^, respectively; thus, the black soliton 
i> = \A^tanh( A /rao'z) exp(— ifit), 



(37) 



c s = vW 



(38) 




(39) 



17 



and velocity y^sin^-. To make a connection to the dark soliton solutions of the NLS 
equation presented in the previous Section, we note that the dark soliton of Eq. fl33|) 
corresponds to a single eigenvalue A = y^sin^. 

Although the system of ZS Eqs. (1391) is linear, its general solution for arbitrary 
initial condition is not available. Thus, various methods have been developed for the 
determination of the spectrum of the ZS problem, such as the so-called quasi-classical 
method [HI IS] (see also Ref. [IS]), the variational approach [ 105] , as well as other 
techniques that can be applied to the case of dark soliton trains [106|,I107] . In any case, 
the generation of single- as well as multiple-dark solitons (see Sec. 13.61 below) can be 
studied in the framework of the 1ST method, and many useful results can be obtained. 
In that regard, first we note that a pair of dark solitons — corresponding to a discrete 
eigenvalue pair in the associated scattering problem — can always be generated by an 
arbitrary small dip on a background of constant density [H] (see also Ref. [15]). This 
means that the generation of dark solitons is a thresholdless process, contrary to the 
case of bright solitons which are created when the number of atoms exceeds a certain 
threshold |108j . In another example, as dark solitons are characterized by a phase jump 
across them, we may assume that they can be generated by an anti-symmetric initial 
wave function profile of the form, 

ip(Zj 0) = y/n^ tanh(az), (40) 

characterized by a background density n and a width a -1 (the ratio ^Jn^/a is assumed 
to be arbitrary). In such a case, the ZS eigenvalue problem f l39|) can be solved 
exactly [H24T8] and the resulting eigenvalues of the discrete spectrum are given by 
Ai = and \ 2 j = —\2j+i = \Jn Q — /i^, where positive fij are defined as fij = y/no — ja, 
j = 1,2, ■ ■ ■ ,N Q , and N is the largest integer such that N < ^/n^/a. These results 
show that for arbitrary ^fn^j a the initial wave function profile of Eq. ( 14*01 will always 
produce a black soliton [cf. Eq. (|37|) ] at z = (corresponding to the first, zero eigenvalue) 
and additional Nq pairs of symmetric gray solitons (corresponding to the even number 
of the secondary, nonzero eigenvalues), propagating to the left and to the right of the 
primary black soliton. Apparently, the total number of eigenvalues and, thus, the total 
number of solitons, is 2N + 1 and depends on the ratio ^Jn^/a. Apart from the above 
example, dark soliton generation was systematically studied in Ref. [15] for a variety 
of initial conditions (such as box-like dark pulses, phase steps, and others). Notice 
that, generally, initial wave function profiles with odd symmetry will produce an odd 
number of dark solitons, while profiles with an even symmetry (as, e.g., in the study of 
Ref. [TJ]) produce pairs of dark solitons; this theoretical prediction was also confirmed 
in experiments with optical dark solitons |109j . Furthermore, the initial phase change 
across the wave function plays a key role in dark soliton formation, while the number 
of dark solitons that are formed can be changed by small variations of the phase. 
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3.3. Integrals of motion and basic properties of dark solitons. 



Let us now proceed by considering the integrals of motion for dark solitons. Taking 
into regard that Eqs. (|271)-(l29l) refer to both the background and the soliton, one may 
follow Refs. [27 t l28 |[TT0|llll] . and renormalize the integrals of motion so as to extract 
the contribution of the background [see Eqs. ( 1321) ]. This way, the renormalized integrals 
of motion become finite and, when calculated for the dark soliton solution fl33|) . provide 
the following results (for k — 0). The number of particles A^s of the dark soliton reads: 



N ds = / (n -\^\ 2 )dz = 2^n- B. (41) 

J — oo 

The momentum Pds of the dark soliton is given by, 



i 



P ds = - I (ipd z il)* - ip*d z ip) dz - n o A0 



/ 00 / 



2 



2 



(42) 



2v{c 2 s -v 2 ) 1/2 + 2c^tan- 1 _ 

where A0 is given by Eq. (135|) and c s = y/n^ is the speed of sound. Furthermore, the 
energy E^ of the dark soliton is given by, 



Erie = — 

ds 2 



\d^\ 2 +{\^\ 2 -n ) 2 ]dz = \{c 2 s -v 2 fl 2 1 (43) 



while the renormalized Lagrangian density takes the form 

c As = l - mr - w) ( i - 1 [\d^\ 2 + {\ip\ 2 - n ) 2 ] . (44) 



2 yr - ^' y [-0)2 y 2 

The renormalized integrals of motion can now be used for a better understanding of 
basic features of dark solitons. To be more specific, one may differentiate the expressions 
f|42|) and f T43|) over the soliton velocity v = A^/n^ to obtain the result, 

^ = (45) 

which shows that the dark soliton effectively behaves like a classical particle, obeying a 
standard equation of classical mechanics. Furthermore, it is also possible to associate 
an effective mass to the dark soliton, according to the equation m ds = dP ds /dv. This 
way, using Eq. ( 1421) . it can readily be found that 

m As = -4^5, (46) 

which shows the dark soliton is characterized by a negative effective mass. The same 
result, but for almost black solitons {B 1) with sufficiently small soliton velocities 
{y 2 <C c 2 ), can also be obtained using Eq. (l4Ul) |112j : in this case, the energy of the dark 
soliton can be approximated as E ds ps (4/3)c^ — 2c s v 2 or, equivalently, 

E ds = E + ^m ds v 2 , (47) 

where E = |c^, and the soliton's effective mass is m ds = — 4 v /n ~. 
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3.4- Small- amplitude approximation: shallow dark solitons as KdV solitons. 

As mentioned above, the case of B 2 = cos 2 1 corresponds to a small-amplitude 
(shallow) dark soliton, which travels with a speed v close to the speed of sound, i.e., 
?; » c s . In this case, it is possible to apply the reductive perturbation method |113] and 
show that, in the small-amplitude limit, the NLS dark soliton can be described by an 
effective KdV equation (see, e.g., Ref. |114j for various applications of the KdV model). 
The basic idea of this, so-called, small- amplitude approximation can be understood in 
terms of the similarity between the KdV soliton and the shallow dark soliton's density 
profile: indeed, the KdV equation for a field u(z,t) expressed as, 

d t u + 6ud z u + d 3 z u = 0, (48) 

possesses a single soliton solution (see, e.g., Ref. [13]): 

u(z, t) = 2K 2 sech 2 [K(z - AkH)} (49) 

(with k being an arbitrary constant), which shares the same functional form with the 
density profile of the shallow dark soliton of the NLS equation [see Eqs. (H9|) and ( 136]) ]. 
The reduction of the cubic NLS equation to the KdV equation was first presented in 
Ref. [9] and later the formal connection between several integrable evolution equations 
was investigated in detail |115j . Importantly, such a connection is still possible even in 
cases of strongly perturbed NLS models, a fact that triggered various studies on dark 
soliton dynamics in the presence of perturbations (see, e.g., Refs. |116till8j for studies in 
the context of optics, as well as the recent review [1 19j and references therein). Generally, 
the advantage of the small-amplitude approximation is that it may predict approximate 
analytical dark soliton solutions in models where exact analytical dark soliton solutions 
are not available, or can only be found in an implicit form [116] . 

Let us now consider a rather general case, and discuss small-amplitude dark solitons 
of the generalized NLS Eq. (fl8l) ; in the absence of the potential (V(z) = 0), this equation 
is expressed in dimensionless form as: 

idtt/> = ~^ + f(n)ip, (50) 

where the units are the same to the ones used for Eq. fT24l) . Then, we use the Madelung 
transformation ip(z,t) = \J n(z, t) exp[i(p(z, t)] (with n = |^| 2 and (p representing the 
BEC density and phase, respectively) to express Eq. fl50|) in the hydrodynamic form: 

d t¥ > + f(n) + \ {d zV f - \n'l 2 d 2 z n'l 2 = 0, (51) 

d t n + d z (nd z ip) = 0. (52) 

The simplest solution of Eqs. f l5T]) - fl5^|) is n = n = \ip \ 2 and = — fit = —fot, where 
/ = f(n ) = f(\ip \ 2 ). Note that in the model of Eq. (JTHJ one has f = ^7=^, for the 
model of Eq. f l20|) . f = ^1 + 2n , and so on. Next, assuming slow spatial and temporal 
variations, we define the slow variables 

Z = e 1/2 (z-ct), T = e 3/2 t } (53) 
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where e is a formal small parameter (0 < e Cl) connected with the soliton amplitude. 
Additionally, we introduce asymptotic expansions for the density and phase: 

n = n + em(Z, T) + e 2 n 2 (Z, T) -\ , (54) 

<P = ~ fot + e 1/2 Vi(Z, T) + &\ % {Z, T) + ■ • • . (55) 

Then, substituting Eqs. (I54l)-(l55l) into Eqs. (I51~])-(l52]). and Taylor expanding the 
nonlinearity function /(n) as f(n) — f + ef Q n x + e 2 [(l/2)/Qn 2 + /on 2 ] + 0( eS ) (where 
/q = g^| n=no ), we obtain a hierarchy of equations. In particular, Eqs. (15T]) - (152"]) lead, 
respectively, at the order 0(e) and 0(e 3 / 2 ), to the following linear system, 

- cd Z (fi + f' ni = 0, n d 2 z (fi - cd z n x = 0. (56) 

The compatibility condition of the above equations is the algebraic equation c 2 = /qUq, 
which shows that the velocity c in Eq. (|5"3l) is equal to the speed of sound, c = c s . 
Additionally, Eqs. (156]) connect the phase <pi and the density ri\ through the equation: 

dz<Pi = —nx. (57) 
To the next order, viz. 0(e 2 ) and C(e 5/2 ), Eqs. and (1521). respectively, yield: 

d T ipi - c s d z p2 + fon 2 + g/o 71 ? + 2 (^^i) 2 - ^ n o ld z n i = °> ( 58 ) 

9 T ni - c s <9 z n2 + <9 Z (rixdzfx) + n d 2 z Lp 2 = 0. (59) 

The compatibility conditions of Eqs. ( 155|) -( 155|) are the algebraic equation c 2 = /qUq, 
along with a KdV equation [see Eq. (PfB]) ] for the unknown density ni: 

2c s d T nx + (3fo + nof^nxdznx - X -d\n x = 0. (60) 

Thus, the density n x of the shallow dark soliton can be expressed as a KdV soliton [see 
Eq. ([4"9~]) ]. In terms of the original time and space variables, n x is expressed as follows: 

*™=- m£*jgf ^ [ ^ K{ *-« )h (61) 

where k is (as before) an arbitrary parameter [assumed to be of order 0(1)], while v is 
the soliton velocity; the latter, is given by 

K 2 

^ = c s -e— , (62) 

and, clearly, v < c s . Apparently, Eq. ( 161]) describes a small- amplitude dip [of order 0(e) 
- see Eq. (154"]) ] on the background density of the condensate, with a phase <pi that can 
be found using Eq. (157]) : in terms of the variables z and t, the result is: 

<fx(z, t) = —- TTTTTT" 777T tanh [ el/2ft; ~ vt)] . (63) 

The above expression shows that the density dip is accompanied by a tanh-shaped 
phase jump. Thus, the wave function characterized by the density nx in Eq. ( 16T]) and 
the phase ipx m Eq. (163]) is an approximate shallow dark soliton solution of the GP 
Eq. (jSip, obeying the effective KdV Eq. (ED}. 
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Notice that the above analysis applies for f{n) = n (i.e., for the cubic NLS model), 
as well as for all forms of the nonlinearity function in Eqs. (|19p -( l2"2"]) . Furthermore, 
variants of the reductive perturbation method have also been applied for the study of 
matter-wave dark solitons in higher- dimensional settings [120)1121] . multi-component 
condensates [122U123] (see also Sec. 16. 1 p and combinations thereof |124j . 

3. 5. On the generation of matter-wave dark solitons 

Matter-wave dark solitons can be created in experiments by means of various methods, 
namely the phase-imprinting, density- engineering, quantum-state engineering (which is a 
combination of phase-imprinting and density engineering), the matter-wave interference 
method and by dragging an obstacle sufficiently fast through a condensate. In connection 
to Sec. 13.21 — and following the historical evolution of the subject — here we will discuss 
the phase-imprinting, density-engineering and quantum-state engineering methods (the 
remaining two methods will be presented in Sees. 16.21 and 16.31 below). 

3.5.1. The phase-imprinting method. The earlier results of Sec. 13.21 as well as more 
recent theoretical studies in the BEC context [125(1126] (see also Ref. |127j ). paved the 
way for the generation of matter-wave dark solitons by means of the phase-imprinting 
method. This technique was used in the earlier [^HISSIHH] - but also in recent [671168] 
- matter-wave dark soliton experiments. The phase-imprinting method involves a 
manipulation of the BEC phase, without changing the BEC density, which can be 
implemented experimentally by illuminating part of the condensate by a short off- 
resonance laser beam (i.e., a laser beam with a frequency far from the relevant atomic 
resonant frequency — see details in the review [128] ). This procedure can be described 
in the framework of Eq. i fPTj) . by considering a time-dependent potential of the form 
V(z;t) oc (p(z)f(t), where f(t) is the laser pulse envelope and <f>(z) is the imprinted 
phase, given by [129J, 



where A</> is the phase gradient, while the width W of the potential edge sets the 
steepness of the phase gradient at z*. Note that since experimentally relevant values 
correspond to a 10-20% absorption width of the phase step, an empirical factor b = 0.45 
is also introduced in Eq. (164")) [129J. 

From a theoretical standpoint, phase- imprinting can be studied (in the absence 
of the trapping potential) in the framework of the 1ST method, upon considering an 
initial wave function of the form ip(z,0) = exp[z</>(z)]; here the imprinted phase <p(z) is 
assumed to increase from left to right and approach constants as z — > ±oo [126] [as, e.g., 
in Eq. ([61)1 ]. The pertinent ZS eigenvalue problem can be solved by mapping Eqs. ( 1391 
to a damped driven pendulum problem. This way, a formula for the number of both 
the even and the odd number of generated dark solitons, traveling in both directions, 
can be derived analytically. 




(64) 
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In some experiments (see, e.g., Ref. [IS]), the generation of the "dominant" dark 
soliton is followed by the generation of a secondary wave packet traveling in the opposite 
direction with a velocity near the speed of sound. This effect can also be understood in 
the framework of 1ST: small perturbations of the dark soliton produce shallow "satellite" 
dark solitons moving with velocities v < c s [T8] . 

3.5.2. The density- engineering method. The density-engineering method involves a 
direct manipulation of the BEC density, without changing the BEC phase, such that 
local reductions of the density are created which eventually evolve into dark solitons. 
This technique was used in the Harvard experiments [I7J,[65]> where a compressed pulse 
of slow light was used to create a defect on the condensate density. This defect induced 
the formation of shock waves that shed dark solitons (or other higher-dimensional 
topological structures, such as vortex rings [65] )■ Notice that the use of a compressed 
pulse of slow light is not really necessary or beneficial in order to create dark solitons 
by means of the density engineering method: in fact, a local reduction of the BEC 
density can also be created by modifying the (harmonic) trapping potential with an 
additional barrier potential, which may be induced by an optical dipole potential or a 
far-detuned laser beam; this barrier can then be switched off non-adiabatically (while the 
harmonic trap is kept on), creating the desired local reduction of the density |129j . This 
technique was employed in a recent experiment [72] , where such a dipole beam was used 
in different setups to induce merging and splitting rubidium condensates; depending 
on the parameters, this process leads to the formation of dark soliton trains, or a high 
density bulge and dispersive shock waves. 

As in the case of phase-imprinting, the density-engineering technique can be studied 
by means of the 1ST method (in the absence of the trapping potential). In fact, earlier 
works [m[T5] (see also Ref. |107j ) have already addressed the problem of dark soliton 
generation induced by initial change of the density: for example, in the case of a box-like 
initial condition, namely i/j(z, 0) = for \z\ > Zq and if)(z, 0) = ^fn{ for \z\ < Zq (with 
Tii < n o)y the ZS spectral problem admits an explicit solution, as it can be solved exactly 
on the intervals \z\ < z and \z\ > z . In particular, it can be shown that there appear 

two discrete eigenvalues Ai,2 = ±2 v /rio [l — < ^ z \{\f n ~^ — a/"T) 2 ] (f° r V^o ~ y/™i ^ v^o) 
and, thus, two small-amplitude dark solitons are generated. 

3.5.3. The quantum-state engineering method. A combination of the phase-imprinting 
and density-engineering methods is also possible, leading to the so-called quantum-state 
engineering technique [129,130j. This method, which involves manipulation of both the 
BEC density and phase, has been used in experiments at JILA [48J and Hamburg [67] 
with a two- component 87 Rb BEC (see Sec. 16. II below): in the one component, a so-called 
"filled" dark soliton was created, with the hole in this component being filled by the 
other component. Depending on the trap geometry, the created filled dark soliton was 
found to be either unstable or stable. Particularly, in the JILA experiment [18] . the dark 
soliton evolved in a quasi-spherical trap (after the filling from the other component was 
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selectively removed) and, due to the onset of the so-called snaking instability, the soliton 
was found to decay into vortex rings (see Sec. 15.11 below). On the other hand, in the 
Hamburg experiment [67], the filled dark soliton in the one component was allowed to 
evolve (in the presence of the other component) in an elongated cigar-shaped trap; this 
way, a so-called dark-bright soliton pair was created (see Sec. 16.11) . which was found to 
be stable, performing slow oscillations in the trap as predicted in theory [131j . 

Notice that a similar two-component engineering technique was also used for the 
creation of vortices [132] . while earlier experimental results from the JILA group could 
be interpreted as a formation of a stack of filled dark solitons in a single BEC [133J. 



3.6. Multiple dark solitons and dark soliton interactions 

3.6.1. The two-soliton state and dark soliton collisions. Apart from the single dark 
soliton solution, the NLS Eq. (|26|) possesses exact analytical multiple dark soliton 
solutions, which can be found by means of the 1ST [12], [21)] (see also Refs. [2T1I22] ). 
Such solutions describe the elastic collision between dark solitons as, in the asymptotic 
limit of t — > ±oo, the multiple-soliton solution can be expressed as a linear superposition 
of individual single-soliton solutions, which remain unaffected by the collision apart from 
a collision-induced phase-shift. To be more specific, let us consider the two-soliton wave 
function ip = tp(z,t), which can be asymptotically expressed as: 

ip ->■ ip(z - ^Jn^Ait,A\,z±) +ip(z - ^fn^A 2 t,A 2 ,z^), t -> +oo, (65) 

ip ->■ ip(z - y/n^Ait,Ai,z^) +ip(z - s /n^A 2 t,A 2 ,z 2 ), t ->• -oo, (66) 

where z 12 denote the position of each individual soliton (in the above expressions, the 
parameters A,- and Bj (j = 1, 2), with A 2 - + B 2 = 1, characterize the velocity and depth 
of the soliton j). Apparently the shape and the parameters of each soliton are preserved, 
while the phase-shift of each soliton is given by: 

■(A 1 -A 2 ) 2 + (B 1 + B 2 ) 2 ^ 



Azi = z7 — z-x — In 

1 1 2Bi 

Az 2 = zt — Zo = In 

2 2 2Bo 



(A, - A 2 ) 2 + (B, - B 2 ] 



(A 1 -A 2 ) 2 + ( J B 1 + J B 2 ) 2 



[(A 1 -A 2 ) 2 + {B 1 -B 



2) 



V2 



(67) 



(6* 



Note that if the soliton velocities are equal, i.e., A\ = —A 2 = A (hence, B\ = B 2 = B), 
then the phase-shift is equal for both solitons and is given by Az = (2B)~ X ln(l+B 2 /A 2 ). 

Equations (!67j) - (!68l) show that the spatial shift of each soliton trajectory is in 
the same direction as the velocity of each individual soliton and, thus, the dark solitons 
always repel each other. Here it should be mentioned, however, that this important result 
(as well as the collision dynamics near the collision point) can better be understood upon 
studying the explicit form of the two-soliton wave function rather than its asymptotic 
limit considered above. To do so, we consider again the case of a two-soliton solution, 
assuming for simplicity that the two solitons are moving with equal velocities (i.e., 
Ai = —A 2 = A). In such a case, the two-soliton wave function is given by 



1>(z,t) = ?pQexp(-itit), (69) 

Lr{Z, t) 
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Figure 2. (color online) The top panel shows the density profile of the two-soliton 
solution in Eq. (|d^|) . The bottom panels show the density profile of two dark solitons 
at their collision point corresponding to z — Zq and t — 0. The density of low-speed 
solitons, v < v c , is characterized by two distinguishable minima (bottom left panel), 
while in the case of high-speed solitons, v > v c , the density exhibits a single minimum 
(bottom right panel); in the critical case, v — v c , the density has a flat single minimum 
(bottom middle panel). 



where 

F(z, t) = 2(n — 2n rnin ) cosh(2noAE>t) 

- 2n A cosh(2^Bz) + i sinh(2n ABt), (70) 
G(z,t) = 2y/n^ cosh(2n ABt) + 2 y /n^cos\i(2^/n^Bz) 1 (71) 

while n m i n = n — n B 2 = n Q A 2 is the minimum density (i.e., the density at the center 
of each soliton). The density profile of the two-soliton solution in Eq. fl69|) is sketched 
in the top panel of Fig. |2j 

To study analytically the interaction and collision between dark solitons, we follow 
the approach of Ref. [71] and find, at first, the trajectory of the soliton coordinate zo as 
a function of time: using the auxiliary equation d z \ip\ 2 = f [where the density \ip\ 2 is 
determined by Eq. f )69|) ]. the following result is obtained: 



cosh^^^o) = ,^^cosh(2n ABt) - 2 1 (72) 

V n min V n o cosh(2n ABt) 

Then, Eq. ( 1721) determines the distance 2zq between the two solitons at the point of 
their closest proximity, i.e., the collision point corresponding to t — 0: 



1 -1 / / n l n ri 



z* = - coslr 1 / _ 2 /™ . (73) 

2Vn - n min V V n min V n o J 

This equation [which holds for n m i n /n = v 2 < 1/4, otherwise Eq. (I731 provides a 
complex (unphysical) value for Zq] shows that Zq = for n mi „/n = A 2 = 1/4. Thus, it 

f Recall that the dark soliton coordinate z is the location of the minimum density (see Fig. [5]). 
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is clear that there exists a critical value of the soliton velocity, namely v c = \^/no = hc s , 
which defines two types of dark solitons, exhibiting different behavior during their 
collision: "low-speed" solitons, with v < v c , which are reflectedby each other, and "high- 
speed" solitons, with v > v c , which are transmitted through each other. In fact, as shown 
in the bottom panels of Fig. [2j the density profile of the low-speed (high-speed) two- 
soliton state exhibits two separate minima (a single non-zero minimum) at the collision 
point, namely u(zq,0) = (u(zq,0) 7^ 0) f . In other words, low-speed solitons are, 
in fact, well-separated solitons, which can always be characterized by two individual 
density minima — even at the collision point — while high-speed solitons completely 
overlap at the collision point. According to the nomenclature of Ref. [134] . the collision 
between slow-speed (high-speed) solitons is called "black collision" ("gray collision"), 
since the dark solitons become black (remain gray) at t = 0. Notice that the case of 
gray collision can effectively be described — in the small-amplitude approximation - 
by the collision dynamics of the KdV equation |134j . 

3.6.2. The repulsive interaction between slow dark solitons. Let us now investigate in 
more detail the case of well-separated solitons, which are always reflected by each other, 
with their interaction resembling the one of hard- sphere-like particles. In particular, we 
consider the limiting case of extremely slow solitons, i.e., no/n m j ra = A 2 \, for which 
the soliton separation is large for every time (i.e., Zq 3> 0); in this case, the second term 
in the right-hand side of Eq. (!72|) is much smaller than the first one and can be ignored. 
This way, the soliton coordinate is expressed as: 

z = — — — cosh^ 1 [A -1 cosh(2n ABt)] . (74) 



The above equation yields the soliton velocities: 
dz Jn^ sinh(2n z/ Bt) 

IT = 1 ( 75 ) 
ai J A- 1 cosh 2 (2n zABt) - 1 

which, in the limit t — > 0, become dzo/dt = 0. Thus, as the dark solitons approach each 
other, their depth (velocity) is increased (decreased), and become black at the collision 
point (t = 0), while remaining at some distance away from each other. Afterwards, the 
dark solitons are reflected by each other and continue their motion in opposite directions, 
with their velocities approaching the asymptotic values dz^/dt = ±y/rioA for t — > ±00 
[see Eq. ( I75|) ]. i.e., the velocity values of each individual soliton. 

Next, differentiating Eq. ( 1741) twice with respect to time, and using Eq. (1721) 
(without the second term which is negligible for well-separated solitons), one may derive 
an equation of motion for the soliton coordinate in the form d 2 zo/dt 2 = —dVi nt (zo)/dzo, 
where the interaction potential Vj nt (zo) is given by: 

Mnt(zb) = \ . , 2 ^°f_ R (76) 
I sinh {2y/n Bzo) 



f In the case of solitons moving with the critical velocity, v — v c — \c s , the two-soliton density exhibits 
a "flat" single zero minimum at the collision point (see bottom-middle panel of Fig. [5]). 
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It is clear that V^nt is a repulsive potential, indicating that the dark solitons repel each 
other. If the separation between the dark solitons is sufficiently large (i.e., 2z ^> 1) 
then the hyperbolic sinh function in Eq. ( 1761) can be approximated by its exponential 
asymptote, and the potential in Eq. (1751) can be simplified as: 

V int (z ) « 2n B 2 expi-A^Bzo). (77) 

The latter expression can also be derived by means of a Lagrangian approach [25J. 
Importantly, although the above result refers to a symmetric two-soliton collision, the 
results of Ref. [TTJ show that it is possible to use the repulsive potential ( ITS"]) in the 
cases of non-symmetric collisions — using an "average depth" of the two solitons - 
and multiple dark solitons — with each soliton interacting with its neighbors (see also 
relevant discussion in Sec. 15.4ft . 

3.6.3. Experiments on multiple dark solitons. Multiple dark solitons were first created 
in a 23 Na BEC in the NIST experiment [46] by the phase-imprinting method (see 
Sec. 13. 5. Tj) . while the interaction and collision between two dark solitons in a 87 Rb BEC 
was first studied in the Hannover experiment of Ref. [19]. Nevertheless, in this early 
experiment the outcome of the collision was not sufficiently clear due to the presence of 
dissipation caused by the interaction of the condensate with the thermal cloud. In the 
more recent Hamburg experiment [68], the phase- imprinting method was also used to 
create two dark solitons in a 87 Rb BEC with slightly different depths. These solitons 
propagated to opposite sides of the condensate, reflected near the edges of the BEC, and 
subsequently underwent a single "gray" collision near the center of the trap. In addition, 
in the recent Heidelberg experiment [69] two dark solitons were created in a 87 Rb BEC 
by the so-called interference method (see Sec. 16.21 below). The solitons observed in 
this experiment, which were "well-separated" ones, propagated to opposite directions, 
reflected and then underwent multiple genuine elastic "black" collisions, from which 
the solitons emerged essentially unscathed. Notice that the experimentally observed 
dynamics of the oscillating and interacting dark-soliton pair of Ref. [69J , as well as the 
one of multiple dark solitons in another Heidelberg experiment [7T], was in a very good 
agreement with theoretical predictions based on the effective particle-like picture for 
dark solitons (see Sec. 14.21 and Sec. 15.41 below) and the interaction potential of Eq. ( 1761) . 

3.6.4- Stationary dark solitons in the trap. At this point, it is relevant to briefly discuss 
the case where multiple dark solitons are considered in a trapped condensate. In this 
case, both the single dark soliton and all other multiple dark soliton states can be 
obtained in a stationary form from the non-interacting (linear) limit of Eq. ( 124"1) . i.e., in 
the absence of the nonlinear term. In this case, Eq. (I24p is reduced to a linear Schrodinger 
equation for a confined single-particle state. For the harmonic potential of Eq. (|9]), 
this Schrodinger equation describes the quantum harmonic oscillator, characterized by 
discrete energies and corresponding localized eigenmodes in the form of Hermite-Gauss 
polynomials [52]. As shown in Refs. [50j[5T], all these eigenmodes exist also in the fully 
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nonlinear problem, and describe an analytical continuation of the above mentioned linear 
modes to a set of nonlinear stationary states. Additionally, analytical and numerical 
results of the recent work [135J suggest that in the case of a harmonic trapping potential 
there are no solutions of the ID GP Eq. (|24j) without a linear counterpart. This actually 
means that interatomic interactions (i.e., the effective mean- field nonlinearity in the GP 
model) transforms all higher-order stationary modes into a sequence of stationary dark 
solitons confined in the harmonic trap [501151]; note that as concerns its structure, this 
chain of, say n, stationary dark solitons shares the same spatial profile with the linear 
eigenmode of quantum number n. From a physical point of view, multiple stationary 
dark soliton states exist due to the fact that the repulsion between dark solitons is 
counter-balanced by the restoring force induced by the trapping potential. 

4. Matter-wave dark solitons in quasi-lD Bose gases 

4-1- General comments. 

We consider again the quasi-lD setup of Eq. ( 1241) . but now incorporating the external 
potential V(z). In this setting, the dynamics of matter- wave dark solitons can be studied 
analytically by means of various perturbation methods, assuming that the trapping 
potential V(z) is smooth and slowly-varying on the soliton scale. This means that in 
the case, e.g., of the conventional harmonic trap [cf. Eq. ()25l)]. the normalized trap 
strength is taken to be Q ~ e, where e 1 is a formal small (perturbation) parameter. 
In such a case, Eq. f )24|) can be expressed as a perturbed NLS equation, namely, 

idtip + \d 2 ^ - |VfV> = R(4>) = V(z)4>. (78) 

Then, according to the perturbation theory for solitons |136j . one may assume that a 
perturbed soliton solution of Eq. (J7SJ) can be expressed in the following general form, 

ip(z,t)=ip s (z,t)+eip r (z,t). (79) 

Here, i[) 8 (z,t) has the functional form of the dark soliton solution (1331) . but with the 
soliton parameters depending on time, and ip r is the radiation — in the form of sound 
waves — emitted by the soliton. Generally, the latter term is strong only for sufficiently 
strong perturbations (see, e.g., Refs. |137[ll38j . as well as Ref. [73] and discussion 
in Sec. I4.4p . Thus, the simplest possible approximation for a study of matter- wave 
dark solitons in a trap corresponds to the so-called adiabatic approximation of the 
perturbation theory for solitons |136j . namely ip(z,t) « tp s (z,t). In any case, the study 
of matter-wave dark solitons in a trap should take into regard that the trap changes the 
boundary conditions for the wave function, and BEC density, namely n — > [instead of 
n — > no in the homogeneous case — see, e.g., Eq. (|36|) ] as z — > ±oo. From a physical 
viewpoint, and based on the particle-like nature of dark solitons (see Sec. I3.3H . one should 
expect that dark solitons could be reflected from the trapping potential; apparently, such 
a mechanism should then result in an oscillatory motion of dark solitons in the trap. 
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There exist many theoretical works devoted to the oscillations of dark solitons in 
trapped BECs. The earlier works on this subject reported that solitons oscillate in 
a condensate confined in a harmonic trap of strength Q, and provided estimates for 
the oscillation frequency. In particular, in Ref. [139j soliton oscillations were observed 
in simulations and a soliton's equation of motion was presented without derivation; 
in the same work, it was stated that the solitons oscillate with frequency Q (rather 
than the correct result which is Vl/\/2 — see below). The same result was derived in 
Ref. [140] . considering the dipole mode of the condensate supporting the dark soliton. 
Other works [141H143] also considered oscillations of dark solitons in trapped BECs. An 
analytical description of the dark-soliton motion, and the correct result for the soliton 
oscillation frequency, fl/v2, were first presented in Ref. [144] by means of a multiple- 
time-scale boundary-layer theory (this approach is commonly used for vortices [53]). 
The same result was obtained in Refs. [11211145] by solving the BdG equations (for 
almost black solitons performing small-amplitude oscillations around the trap center 
- see Sec. 14.31 below), using a time- independent version of the boundary- layer theory. 
Furthermore, in Ref. |112j a kinetic-equation approach was used to describe dissipative 
dynamics of the dark soliton due to the interaction of the BEC with the thermal cloud. 

Matter-wave dark soliton dynamics in trapped BECs was also analyzed in other 
works by means of different techniques that were originally developed for optical dark 
solitons [31]. In particular, in Ref. [146] the problem was analyzed by means of the 
adiabatic perturbation theory for dark solitons devised in Ref. [28], in Ref. [147] by 
means of the small-amplitude approximation (see Sec. I3.4f) . while in Ref. [148] by means 
of the perturbation theory of Ref. [29]. Later, in Refs. [149, 150J the so-called "Landau 
dynamics" approach was developed, based on the use of the renormalized soliton energy 
[cf. Eq. (|43p ], along with a local density approximation. Models relevant to the dynamics 
of matter-wave dark solitons in ID strongly-interacting Bose gases, were also considered 
and analyzed by means of the small-amplitude approximation [15111152] (see also work for 
dark solitons in this setting in Refs. [153H157] ). In other works, a Lagrangian approach 
for matter- wave dark solitons was presented [158] (see also Ref. [159] ). and an asymptotic 
multi-scale perturbation method was used to describe dark soliton oscillations and the 
inhomogeneity-induced emission of radiation [160J. Recently, the motion of dark solitons 
was rigorously analyzed in Ref. [161] (where a wider class of traps was considered), while 
in Ref. [162] the same problem was studied in the framework of a generalized NLS model. 

Finally, as far as experiments are concerned, the oscillations of dark solitons were 
only recently observed in the Hamburg [671(68] and Heidelberg [69|I7T] experiments. In 
these works, the experimentally determined soliton oscillation frequencies were found to 
deviate from the theoretically predicted value Cl/y/2. This deviation was explained 
in Refs. J67J[68] by the anharmonicity of the trap, while in Refs. |69. 71] by the 
dimensionality of the system and the soliton interactions (see also Sec. 15.41 below). 
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4-2. Adiabatic dynamics of matter-wave dark solitons 

4-2.1. The perturbed NLS equation. The adiabatic dynamics of dark matter- wave 
solitons may be studied analytically by means of the Hamiltonian [27J 28J or the 
Lagrangian [25] approach of the perturbation theory for dark solitons, which were 
originally developed for the case of a constant background. These approaches were 
later modified (see Ref. |146] for the Hamiltonian approach and Refs. |158[I159] for the 
Lagrangian approach) to take into regard that, in the context of BECs, the background 
is inhomogeneous due to the presence of the external potential. The basic steps of these 
perturbation methods are: (a) determine the background wave function carrying the 
dark soliton, (b) derive from Eq. (ITS]) a perturbed NLS equation for the dark soliton 
wave function, and (c) determine the evolution of the dark soliton parameters by means 
of the renormalized Hamiltonian [cf. Eq. ( H3l) ] or the renormalized Lagrangian [cf. 
Eq. (JED] of the dark soliton. Here, we will present the first two steps of the above 
approach and, in the following two subsections, we will describe the adiabatic soliton 
dynamics in the framework of the Hamiltonian and Lagrangian approaches. 

We consider again Eq. (ITS"]) and seek the background wave function in the form, 

ip(z, t) = $0) exp(-i/it + i6 ), (80) 

where // is the normalized chemical potential, 6 Q is an arbitrary phase, while the unknown 
real function $(z) satisfies the following equation, 

2 ~cb 2 

Then, we seek for a dark soliton solution of Eq. (ITS"]) on top of the inhomogeneous 
background satisfying Eq. (18T1) . namely, ip = $(z) exp(— z/it + i9 )tp s (z,t), where the 
unknown wave function ip s (z, t) represents a dark soliton. This way, employing Eq. (1ST]) , 
the following evolution equation for the dark soliton wave function is readily obtained: 

iW. + X -dl^ s - $ 2 (|^| 2 - l)ip s = ~ hi($)<9^. (82) 

It is clear that if the trapping potential V(z) is smooth and slowly- varying on the soliton 
scale, then the right- hand- side, and also part of the nonlinear terms of Eq. (1821) . can be 
treated as a perturbation. To obtain this perturbation in an explicit form, we use the 
TF approximation to express the background wave function as &(z) = ^/l — V{z) [see 
Eq. (TlOj) for g = 1 and /x = 1 f ] and approximate the logarithmic derivative of $ as 

_ ± M ~^h +V + V 2 ). (83) 
dz 2 dz K ' K ' 

This way, Eq. (|82|) leads to the following perturbed NLS equation, 

^ + \^-m 2 -i)^ = Q^ s \ (84) 

where the perturbation Q(ip s ) is given by: 

Q(^) = (1 - |^| 2 ) 1> 9 V + \d z ^(l + V + V 2 ). (85) 

f It can easily be shown that the main result of the analysis [cf. Eq. (|90|) ] can be generalized for every 
value of n such that the system is in the TF-1D regime. 



^+ o^TT-^ 3 = V{z)$. (81) 
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4-2.2. Hamiltonian approach of the perturbation theory. First we note that in the 
absence of the perturbation fl85|) . Eq. f[Mj) has a dark soliton solution of the form: 

^ s (z, t) = cos0tanh£ + isin0, (86) 

where ( = cos0[x — (sin0)t] [see Eq. ([33]) ] . Then, considering an adiabatic evolution 
of the dark soliton, we assume that in the presence of the perturbation the dark soliton 
parameters become slowly- varying unknown functions of t [27 1 128 | [T46] . Thus, the 
soliton phase angle becomes — > <fi(t) and, as a result, the soliton coordinate becomes 
C ~~ y C(t) = cos 4>(t) [z — z (t)}. In the latter expression, the dark soliton center z (t) is 
connected to the soliton phase angle through the following equation: 

^ -*.«*). (87) 

The evolution of the soliton phase angle can be found by means of the evolution of the 
renormalized soliton energy. In particular, employing Eq. f [4"3"l) (for /i = 1), it is readily 
found that dE^s/dt = — 4 cos 2 (ft sin (j)(d(fi / dt) . On the other hand, using Eq. ([84"[) and its 
complex conjugate, it can be found that the evolution of the renormalized soliton energy 
is given by dE ds /dt = — f_°° [Q(ip s )d t ip* + Q*(ip s )dt4> s } dz. Then, the above expressions 
for dE&s/dt yield the evolution of <f), namely 

' " Q(i> s )d t r s dz 



, . -Re / QMd t r s dz . (88) 
at 2 cos z (p sin |_ J-oo 

Next, we Taylor expand the potential V(z) around the soliton center zq, and assume 
that the dark soliton is moving in the vicinity of the trap center, i.e., \i = 1 ^> V, which 
means that the last two terms in the right-hand side of Eq. (JHSJ) can be neglected. This 
way, one may further simplify the expression for the perturbation in Eq. ([85]) which, 
when inserted into Eq. ( [881) . yield the following result: 



^ = -2 COS< W (89) 
To this end, combining Eq. ([89]) with Eq. f[H7j) . we obtain the following equation of 
motion for nearly stationary (black) solitons with cos0 ~ 1, 

^ = -i^ (90) 

The above result indicates that the dark soliton center can be regarded as a Newtonian 
particle: Eq. ( 190]) has the form of a Newtonian equation of motion of a classical particle, 
of an effective mass M e fj = 2, in the presence of the external potential V. In the case of 
the harmonic potential [cf. Eq. ( 125]) ]. Eq. ( 190]) becomes the equation of motion of the 
classical linear harmonic oscillator, d 2 z Q /dt 2 = — (1/2)Q 2 zq, and shows that the dark 
soliton oscillates with frequency 

Uosc = -J=, (91) 

or, in physical units, with co z / \/2. An example of an oscillating matter- wave dark soliton 
is shown in Fig. 
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Figure 3. (Color online) Contour plot showing the evolution of the density of a 
harmonically confined BEC, as obtained by direct numerical integration of the GP Eq. 
([75)) . The initial condition is ip = [fi - (1 / 2)9? z 2 ] tanh(z - z (0)), i.e., a TF cloud, 
characterized by a chemical potential fi, and carrying a dark soliton initially placed at 
z = z (0). The parameter values are ^ = 1, fi = 0.025 and z (0) = 4. The dark soliton 
oscillates with frequency oj osc — 17/ a/2 ~ 0.018. The dotted line across the soliton 
trajectory corresponds to the analytical prediction of Eq. (|90|) . 



At this point, it is relevant to follow the considerations of Ref. |112j (see also 
[14411145] ) and estimate the energy E ds of this almost dark soliton in the trap. Taking 
into regard that in the case of a homogeneous BEC this energy is given by Eq. f[47|) . 
one may use a local density approximation and use in Eq. ( [43]) the local speed of sound, 
c(z) = \Jn a {z) [79] (here, n (z) is the density of the ground state of the BEC), rather 
than the constant value c s = y/n^ [cf. Eq. (I38j) ] . Then, in the TF limit, the density is 
expressed as n (z) — ji — \Q 2 z 2 = c 2 — ^ 2 z 2 and, thus, one may follow the lines used 
for the derivation of Eq. f[471) (for sufficiently slow solitons and weak trap strengths) and 
obtain the result: 

£ds = E + ^m ds v 2 + -m dB Q 2 z 2 , (92) 

where E = and m ds = — 4y / n^ as in Eq. f[4T|) . The above equation shows that the 
incorporation of the harmonic trap results in a decrease of the energy of the dark soliton 
by the potential energy term i|md s |^ 2 £ 2 . Moreover, the ratio of the soliton mass over 
this potential energy is given by (Q 2 z 2 /4) _1 , which is exactly two times the ratio of the 
atomic mass (which is equal to m = 1 in our units) over the external potential, namely 
(Q 2 z 2 /2) _1 . This is another interpretation of the result that the effective mass of the 
dark soliton center is M e s = 2. 

4-2.3. Lagrangian approach for matter-wave dark solitons. The perturbed NLS Eq. 
(184|) . with the perturbation of Eq. (1851) . can also be treated by means of a variational 
approach as discussed in the beginning of Sec. 14.21 First, we assume that the solution 
of Eq. ([HID is expressed as [see Eqs. ([33]) and ([341 ]: 

ip s (z,t) = EtanhC + iA (93) 
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Here, A and B are unknown slowly-varying functions of time (with A 2 + B 2 = 1) 
representing, respectively, the velocity and amplitude of the dark soliton (which become 
time-dependent due to the presence of the perturbation), while ( = B(t) [z — Za(t)], 
where z (t) is the dark soliton center. Note that in the unperturbed case, dz /dt = A, 
but in the perturbed case under consideration, this simple relationship may not be valid 
(see below). Next, the evolution of the unknown soliton parameters ctj(t) (which is a 
generic name for z (t) an d A(t)) are obtained via the Euler-Lagrange equations [25|ll58j : 

dL ds _ d_f dL ds 
dt 



dotj 



da. 



2Re 



Q*(^)^d: 

oo 



(94) 



where otj 



dacj/dt and L ds = f_°° dzC ds {ip s } represents the averaged Lagrangian of 
the dark soliton of the unperturbed NLS equation (namely for Q(ip s ) = 0), with the 
Lagrangian density £ds being given by Eq. f|44|) (for no = 1). The averaged Lagrangian 
can readily be obtained by substituting the ansatz ( l93l) into Eq. ( |44l) : 



L 



ds 



dzQ 



-i 



- A -B\ 
3 



(95) 



B 

— AB + tan ^ \ — 

dt \ A, 

Therefore, substituting Eqs. ( 1951) and (1851) into Eq. ( 1941) . it is straightforward to 
derive evolution equations for the soliton parameters. For completeness, we will follow 
Ref. |158j and present the final result taking also into account the last two terms in the 
right-hand side of Eq. (|85|) — which were omitted in the previous subsection — so as 
to describe the motion of shallower solitons as well. This way, and employing a Taylor 
expansion of the potential around the soliton center (as in the previous subsection), we 
obtain the following evolution equations for zo(^) and A(t): 

r 2\ / qtA 2 



dzo 
~dT 



A 



1 



-V(z ) 



A 
4B 2 



7T 
"9 



dV 

dz 



[1 - 2V{zq)] 
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~dt 
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dV 

dz 
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Equations f l96|) -f l97|) describe the dark soliton dynamics in the trap, in both cases of 
nearly black solitons (A m or B « 1) and gray ones (with arbitrary A or B). In the 
former case, and neglecting the higher-order corrections arising from the inclusion of 
the last two terms in the right-hand side of Eq. (185|) . the result of Eq. (|90|) is recovered: 
nearly black solitons oscillate near the trap center with the characteristic frequency 
given in Eq. fl9T]) . On the other hand, numerical simulations in Ref. |158] have shown 
that the full system of Eqs. ( I96"l) -(l97l) predicts that shallow solitons oscillate in the trap 
with the same characteristic oscillation frequency. Therefore, there is a clear indication 
that the oscillation frequency of Eq. (!9~TT) does not depend on the dark soliton amplitude. 
This result is rigorously proved by means of the Landau dynamics approach that will 
be discussed below. 
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4-2.4- Landau dynamics of dark solitons. The oscillations of dark solitons of arbitrary 
amplitudes in a trap can also be studied by means of the so-called Landau dynamics 
approach devised in Refs. [144,150]. This approach, which further highlights the particle- 
like nature of the matter-wave dark solitons, relies on a clear physical picture: when 
a dark soliton moves in a weakly inhomogeneous background, its local energy stays 
constant. Hence, one may employ the local density approximation, and rewrite the 
energy conservation law of Eq. (j4"3"j) as c 2 (z ) — v 2 = (3E ds /4:) 2 ^ 3 , where c 2 (z ) is the 
local speed of sound evaluated at the dark soliton center zq. Then, in the TF limit, one 
has c 2 (z ) = c 2 — \Q 2 z 2 (as before), and taking into regard that the soliton velocity is 
v = dzo/dt, the following equation for the energy of the dark soliton is readily obtained: 



where E ds = c 2 — (3i?ds/4) 2 / 3 and the effective mass of the dark soliton center is again 
found to be M cS = 2. It is readily observed that Eq. f )98|) can be reduced to Eq. f )90|) and, 
thus, it leads to the oscillation frequency of Eq. f l9~Tj) . Nevertheless, the result obtained 
in the framework of the Landau dynamics approach is more general, as it actually refers 
to dark solitons of arbitrary amplitudes. Moreover, the same approach can be used 
also in the case of more general models, including, e.g., the cases of non-harmonic traps 
and/or more general nonlinearity models, such as the physically relevant ones described 
by Eqs. ( fl~9l -(l22l) |150j . Nevertheless, it should be noted that in such more general cases 
the problem can be treated analytically for almost black solitons (v <C c), performing 
small-amplitude oscillations. In this case, the conservation law E ds (c(zo), v ) = Eq can 
be Taylor expanded around z = and v = 0, leading to expressions for the soliton's 
effective soliton mass and oscillation frequency |150j . 

4-2.5. The small- amplitude approximation. Next, we discuss the adiabatic dynamics of 
small- amplitude dark solitons in trapped ID Bose gases. In this case, one may formally 
reduce the more general GP model of Eq. (150]) (including the potential term Vip) to 
a KdV equation with variable coefficients — see Ref. [Sj for details and Ref. |163] 
for applications of this KdV model. The main result of such an analysis is that the 
density and the phase of the approximate shallow dark soliton solution of Eq. ( !50l) 
have, to the leading-order of approximation, the functional form of their counterparts in 
Eqs. ( )6T|) and ( 163]) . but with the soliton parameter n depending on a slow variable, say 
Z (see earlier work for a calculation of k(Z) in Refs. [164, 165J). This way, approximate 
analytical shallow soliton solutions have been found in various works |f 34[ll47lll51U152] 
for different forms of the nonlinearity function f(n). Nevertheless, there are some subtle 
issues concerning the validity of this approximation, as discussed in Refs. [148[I150[I160| . 
which is, strictly speaking, valid away from the turning points (where the soliton velocity 
vanishes). On the other hand, numerical results (see, e.g., Ref. |152j ). illustrate that 
the range of validity of the above results is, in fact, wider than what may be expected 
based on the limitations of this approach. 




(98) 
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The small-amplitude approximation, along with a local- density approximation, has 
also been used to estimate the shallow soliton's oscillation frequency: the shallow 
soliton's velocity v, which in the homogeneous problem was found to be close to the 
speed of sound, namely v ~ c s = a/7o^o [see Eq. fl62l ], can be approximated in the 
inhomogeneous system as follows: 

^c s (Z) = vS(Ij. (99) 

In some cases, Eq. (|99|) can be used for the derivation of physically relevant results. 
For example, following Ref. |15U] . we assume that f(n) = n a , where a = 1 for weakly- 
interacting BECs, or a = 2 for strongly-interacting Tonks gases |101j . Then, in the TF 
limit, n {Z) = \ji - V{Z)] l l a , and Eq. (EH} is reduced to the form dZ / y/fj, - V{Z) = 
y/adt. The latter is integrated and yields [for V(Z) = (1/2)Q 2 Z 2 ] the soliton trajectory: 



Z = LBm[(ny/a/2)t], (100) 

where L = y/2jl/VL is the length of the TF cloud. Equation (jlOOp predicts that the 
shallow dark soliton will perform oscillations approximately in the entire spatial region 
occupied by the gas, with an oscillation frequency which takes the following values (in 
physical units): for a = 1, i.e., for quasi-lD BECs described by the cubic GP equation, 
Cow = w z /V2, while for a = 2, i.e., for the Tonks gas described by a quintic NLS 
equation, u osc = u z . Note that the latter result was first obtained via a many-body 
calculation |166j . and later was derived by means of the KdV approximation |151j . 

4-3. Bogoliubov-de Gennes analysis of stationary dark solitons 

4-3.1. The single- soliton state. The above result, namely the fact that almost black 
solitons perform oscillations around the trap center with the frequency given in Eq. (l9ip. 
can also be derived by means of a BdG analysis as was first demonstrated in Ref. |145] 
(see also results in Refs. |167H169j ). Such an analysis can be done in the TF limit for the 
stationary dark soliton state, namely the black soliton ip [see Eq. (|37|) ] located at the 
trap center, i.e., z$ = 0, which is actually the first excited state of the condensate. Then, 
following the discussion in Sec. 12.31 the excitation spectrum can be found as follows: 
using the ansatz ip(z,t) = [ipo{ z ) + u(z)e~ luJt + v*(z)e tujt ] e~ iflt [where u — oj r + iu>i is a 
(generally complex) eigenfrequency and (u, v ) are perturbation eigenmodes], we derive 
from Eq. ff24|) the following BdG equations: 

[H - n + Vq]« + = LOU, (101) 

[H - /i + ^l\v + i;* 2 u = -uv, (102) 

where H = —(l/2)d 2 + (l/2)Q 2 z 2 is the single particle operator. A typical example 
showing the initial configuration, i.e., the condensate and the stationary dark soliton 
(which can be found, e.g., by a Newton- Raphson method), as well as the corresponding 
spectral plane (co r ,Ui), are shown in Fig. HI 

The BdG analysis reveals that all the eigenfrequencies of the spectrum are real, 
which indicates that the stationary dark soliton is dynamically stable. The four 
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Figure 4. (Color online) Left panel: The density of a condensate carrying a stationary 
(black) soliton located at z = 0. The normalized chemical potential is /i = 1. 
The (green) dashed line shows the trapping potential with normalized trap strength 
ft = 0.025. Right panel: The lowest characteristic eigenfrequencies of the Bogoliubov 
excitation spectrum: The eigenfrequency located at the origin corresponds to the 
Goldstone mode, the one at = 0.025 to the Kohn mode, and the one at y/3n to 
the quadrupole mode. Finally, there exist a anomalous mode with ua — £l/y/2. 

smallest magnitude eigenfrequency pairs f and their corresponding eigenmodes have 
the following physical significance (see, e.g., Ref. [42J). First, there exists a zero 
eigenfrequency, located at the origin of the spectral plane (u r ,Ui), which reflects the 
phase invariance of the ID GP equation. The respective eigenfunction is the so-called 
Goldstone mode and does not result in any physical excitation (oscillation) of the 
system. The solutions with eigenfrequencies u r = ±f2 correspond to the so-called dipole 
mode (or Kohn mode), which is relevant to the motion of the center of mass of the 
system; note that as the system is harmonically confined, the center of mass oscillates 
with the frequency of the harmonic trap |170j . The solutions with eigenfrequencies 
having the next larger magnitude eigenfrequencies correspond to the quadrupole mode, 
with the location of the eigenfrequencies at u r = ±v3$7, being particular to the one- 
dimensionality of the system [88J. Note that the excitation of the quadrupole mode 
(induced, e.g., by a time-modulation of the trap strength) results in a breathing behavior 
of the condensate, with its width oscillating with the above-mentioned frequency. 

Of particular interest are the solutions with eigenfrequencies u r = oja = n/y/2, 
which correspond to the so-called anomalous mode. This mode appears in the 
Bogoliubov analysis only when topological excitations of the condensate are involved, 
namely dark solitons or vortices [74]. A characteristic property of the anomalous mode 
is that the integral of the norm x energy product, /(|w| 2 — \v \ 2 )udz (in our units), is 
negative rather than positive as is the case for all the positive frequency modes associated 

f Recall that due to the Hamiltonian nature of the system, the eigenfrequencies ±u; r correspond to the 
same physical oscillation. 
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Figure 5. (Color online) Left panel: The eigenfunctions ua (solid line) and va (dotted 
line) of the anomalous mode. Right panel: The solid (black) line shows the density 
of the stationary dark soliton in a region around the trap center. The dashed (blue) 
line shows the density of the dark soliton when excited by the anomalous mode. The 
parameter values are the same to the ones of Fig. |U 

with the ground state of the system [42J. Note that, from a mathematical viewpoint, 
the anomalous mode possesses a topological property of the so-called negative Krein 
signature |171j . namely K = sign{J(|w| 2 — |f| 2 )o;(i^} < (for positive eigenfrequencies 
u>). Practically, this means that the anomalous mode becomes structurally unstable 
(i.e., it becomes complex) upon collision with other eigenvalues, as is the case when 
dissipation is present |172j . In our case, finite temperature automatically implies the 
presence of dissipation which, as discussed in more detail in Sec. 16.51 below, may be 
described - - in the simplest possible approach - - by including a phenomenological 
temperature-induced damping term in the GP model. 

In order to better clarify the above, and discuss in more detail the stability of 
the excitation corresponding to the anomalous mode, namely of the dark soliton, we 
note the following. At temperatures T — > (as is the case under consideration), the 
negative energy of the dark soliton does not imply any instability (e.g., a decay process) 
and, thus, the soliton is dynamically stable. Nevertheless, at finite temperatures, i.e., 
in the presence of a thermal cloud, the above mentioned properties of the anomalous 
mode indicate that the soliton will become unstable: in this case, the presence of the 
temperature-dependent damping results in the decay of the soliton (see discussion in 
Refs. [112(1145] as well as in Sec. 16.51 below). From a physical point of view, the decay 
mechanism resembles the one of the low-energy excitations of trapped BECs [173j and 
originates from the scattering of thermal particles on the dark soliton. Thus, according 
to these arguments, matter-wave dark solitons can be regarded as thermodynamically 
unstable excitations as, in the presence of the temperature-induced dissipation, the 
system will be driven towards configurations with lower energy; in other words, the 
dark solitons will decay to the ground state. This scenario is also often referred to as 
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Figure 6. (Color online) Left panel: The density of a condensate carrying a stationary 
two-dark-soliton state (parameter values are fi = 1 and fi = 0.025, as in Fig. @|. 
The solitons are located at z = ±2.3. Right panel: The lowest characteristic 
eigenfrequencies of the Bogoliubov excitation spectrum: The eigenfrequency located 
at the origin corresponds to the Goldstone mode, the one at fi = 0.025 to the Kohn 
mode, and the one at \/3fi = 0.043 to the quadrupole mode. Finally, there exist two 
anomalous modes with eigenfrequencies U\ = 0.0179 and w 2 = 0.0566. 

energetic instability [174]. 

As mentioned above, the eigenfrequency of the anomalous mode is equal to the 
oscillation frequency of a dark soliton around the center of the harmonic trap in the TF 
limit. On the other hand, the eigenfunctions ua and va of the anomalous mode, shown 
in Fig. EJ are localized within the notch of the dark soliton |145[I169] . and their sum 
can be approximated as ua + va oc sech 2 (y/n^z). Notice that in the case of a uniform 
condensate (i.e., in the absence of the trap), there exists a translational mode with zero 
frequency which has the same functional form, namely d z ipo, due to the translational 
invariance of the dark soliton solution. When the trap is present, however, this symmetry 
is broken, which suggests that the anomalous mode can be regarded as the "ghost" of 
the broken translational invariance of the dark soliton solution. We finally mention that 
the direct connection of the anomalous mode to the oscillation of the dark soliton can be 
better explained by the fact that an excitation of the stationary black soliton ip by the 
anomalous mode results in a displacement of the soliton from the trap center. In other 
words, the GP Eq. ( |24|) with the initial condition i/;(z; t = 0) = ijJo{z) + ua{z) + v A (z), 
will naturally lead to dark soliton oscillations studied in Sec. 14.21 

4-3.2. The multi- soliton state. The BdG analysis can also be performed in the more 
general case of multi-soliton states P-75J, which may be found in a stationary form 
(as explained in Sec. 13 .6 j) . In this case, starting from the non-interacting limit, it can 
be found that the Bogoliubov spectrum of the n-th excited state consists of one zero 
eigenvalue (corresponding to the Goldstone mode), n double eigenvalues (accounted for 
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Figure 7. (Color online) Spatio-temporal contour plot of the condensate density, for 
parameter values fj, = 1 and fi = 0.025 (as in Fig. [BJ. Top panel: The dark solitons, 
initially placed at z± — and Z2 — 5, oscillate in-phase with a frequency u; osc = 0.018 ~ 
u)\ = 0.0179 fa f2/y2 = 0.0176. Bottom panel: The dark solitons, initially placed at 
2 = ±3, oscillate out-of-phase with a frequency w osc = 0.057 ~ u>2 = 0.0566. Here, 
oji.2 are the eigenfrequencies of the first and second anomalous mode, respectively. 

by the presence of the harmonic trap), and infinitely many simple ones. Then, in the 
nonlinear regime, one of the eigenvalues of each double pair becomes an anomalous 
mode of the system (characterized by a negative valued integral of the norm x energy 
product) and, thus, the number of anomalous modes in the excitation spectrum equals 
to the number of dark solitons |175j . This is in agreement with the fact that the number 
of eigenvalues with negative Krein signature equals to the number of the nodes of the 
stationary state [176] . Notice that in the framework of the ID GP Eq. (|24|) — i.e., in 
the TF-1D regime — the first anomalous mode coincides with the oscillation frequency 
co osc = Cl/y/2 of the single dark soliton. An example of a condensate with a stationary 
two-dark soliton state, as well as the pertinent spectral plane, are shown in Fig. [6j 

The physical significance of the n-anomalous modes has been discussed in Refs. [HJ 
175J: for example, in the case of a two-dark soliton state, the smallest of the anomalous 
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modes corresponds to an in-phase oscillation (i.e., when the two dark solitons oscillate 
together without changing their relative spatial separation), the largest anomalous mode 
corresponds to an out-of-phase oscillation (i.e., when the two dark solitons move to 
opposite directions with the same velocity and undergo head-on collision), and so on. 
An example of the in-phase and the out-of-phase oscillation of the two-dark soliton 
state shown in the left panel of Fig. [6] (when the solitons are properly displaced from 
their equilibrium positions) is illustrated in Fig. Here it should be noticed that since 
the starting point of the above considerations is the non- interacting limit, a similar 
analysis can also be performed in the case of other mean-field models with non-cubic 
nonlinearity, as the ones describing cigar-shaped BECs in the dimensionality crossover 
regime from 3D to ID (see Sec. 12.41 and Sec. 15.41 below). 

Finally, as concerns the stability of nonlinear modes (see relevant investigations in 
Refs. [5Tll711ll~7~~ 1178] ). all higher-order nonlinear modes are unstable near the non- 
interacting limit [71 ] 1178] (but can be stabilized by using anharmonic traps |178j ); 
nevertheless, the instability ceases to exist sufficiently deep inside the nonlinear regime 
(i.e., for sufficiently large BECs, with large N) [71] . 

4-4- Radiation effects: inhomogeneity-induced sound emission by the soliton. 

As indicated in Sec. 14.11 the dark soliton experiences a background density gradient in 
the presence of the perturbation R{jp) in Eq. ( ITS]) and, thus, it continuously emits energy 
in the form of sound waves. Here, we will study this effect in more detail, considering 
the case of R(ip) = V(z)ip with V(z) being a harmonic potential. A particular feature of 
this setup, is that the emitted sound energy remains confined within the spatial region 
of the trap and, hence, the soliton continuously re-interacts with the emitted sound 
waves |137j ; in fact, this process is such that, on average, the dark soliton reabsorbs 
the radiation it emits. Thus, in the case of harmonic traps, an investigation of the 
inhomogeneity-induced sound emission, as well as an estimation of the rate of emission 
of energy, is relevant for short timescales, i.e., < t < T osc = 27i/u osc [where u osc is given 
by Eq. (}91~) . with O < 1], On the other hand, if dark solitons evolve in the presence 
of non-harmonic potentials (such as localized barriers pT3l .ll38j.n~46l.ll59l.il T9j . disordered 
potentials [180U181 j . anharmonic traps [182], or other "properly designed" potentials - 
see below), the problem may be easier - at least in terms of a numerical investigation: 
in fact, as is explained below, it is possible to consider suitable setups that either damp 
off the emitted sound density or cause the emitted sound to dephase. 

Various such setups were proposed and analyzed in the past; the most prominent 
example is, perhaps, a tight inner "dimple" trap, confining a dark soliton, located within 
a much weaker outer harmonic potential (such a configuration can be realized by focusing 
an off- resonant laser beam within a harmonic trap) |137j . In this case, if the depth of the 
dimple trap is sufficiently small, the sound waves can escape (to the outer trap), while 
the soliton can remain confined in this region. In this limit, sound energy is damped off 
for short enough timescales, until it bounces off the weaker outer trap and thus becomes 
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forced to re-interact with the soliton. An alternative setup considered in Refs. |183]I184] 
consists of a harmonic trap perturbed by an optical lattice potential (see Sec. I6.4p . In 
this case, the optical lattice can confine a soliton within a single or a few lattice sites, 
with the sound (again for short enough times) escaping to neighboring sites. Although 
in this case the sound re-interacts with the soliton on faster timescales than in the case 
of the dimple trap mentioned above, the presence of the lattice dephases the emitted 
sound waves, and hence accelerates the soliton decay. 

The radiation-induced dissipation of matter-wave dark solitons in harmonic traps 
was studied analytically in Ref. |160] by means of an asymptotic multi-scale expansion 
method. Particularly, assuming that R(ip) = e 2 z 2 ip (with e being a formal small 
parameter defined by the aspect ratio Q), the following results were obtained. In 
the limit e — > 0, the dark soliton evolves adiabatically so that the dark soliton center 
Za(t) = vt + z Q — >■ s(T)/e, i.e., it becomes a function of the slow timescale T = et, while 
the soliton velocity is given by v(T) — s = ds/alT. The adiabatic dynamics is followed 
by generation of sound waves, which can be taken into regard as per Eq. ( 1791) . In fact, in 
the decomposition of the wave function ip into an inner and an outer asymptotic scale, 
the leading-order radiative effects are taken into account when the complex phase 9 Q [see 
Eq. (180]) ] depends also on T = et, i.e., O = 9(T), and the first-order corrections to the 
dark soliton (|33|) grow linearly in z. Neglecting reflections from the trap, the extended 
dynamical equation for the position s(T) of the dark soliton (133]) takes the form: 



s + s = —= - = + Q(e 2 ). (103) 

2y/(l ~ S 2 )Vl " S 2 - S 2 

The left-hand-side of Eq. (11031) represents the leading-order adiabatic dynamics of the 
dark soliton [see also Eq. ( 190]) ] oscillating on the ground state of the trap, namely 
a harmonic oscillator with the obvious solution s(T) = socos(T + So) (with s and Sq 
being arbitrary parameters representing the initial position and phase of the soliton). As 
long as s 2 + s 2 < 1, the adiabatic dynamics approximation remains valid for large values 
of the position s(0) and speed s(0) of the dark soliton. In other words, dark solitons 
oscillate in the trap with a uniform oscillation frequency for solitons of all amplitudes 
and velocities, in agreement with the prediction of the Landau dynamics approach (see 
Sec. I4.2p . Apparently, in the limiting case of s 2 + s 2 — > 1 (i.e., for extremely shallow 
dark solitons) Eq. (I103f) is not applicable. 

Next, letting E = \ (s 2 + s 2 ) be the energy of the harmonic oscillator, one may 
employ Eq. (I103P to calculate the rate of change of E due to the leading-order radiative 
effects appearing in the right-hand side of Eq. ( I103p . The result is: 

E = -. + 0(e 2 ) > 0, (104) 

2^(1 - s 2 )Wl - s 2 -s 2 

and shows that due to the energy pumping (11041) . the amplitude of the harmonic 
oscillator increases in time. In the limit s 2 + s 2 — > 0, Eqs. (I103p and (I104p can be 
simplified. First, the energy of the dark soliton oscillations accelerates by the squared 
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law E = es 2 /2, which was confirmed in numerical simulations in the setup of Ref. |137j . 



Second, the nonlinear equation HI 03[) is linearized as follows: 




Equation (I105P includes an anti- damping term accounting for the emission of radiation, 
indicating that the center point (0, 0) becomes an unstable spiral point on the plane 
(s, s). Apparently, the leading-order solution reads s(T) = s e eT//4 cos(T + <5 ); thus, the 
amplitude of oscillations of a dark soliton increases while its own amplitude decreases. 

The above results of the asymptotic analysis were confirmed by the numerical 
simulations of Ref. |160j . but also by numerical findings reported in other works: 
the radiation-induced effects were also observed for dark solitons oscillating between 
two Gaussian humps [138J, or for dark solitons that are parametrically driven by a 
pair of two periodically-modulated Gaussian barriers, oscillating in anti-phase at a 
frequency close to the soliton frequency |185j . It is worth noticing that the mechanism 
proposed in Ref. |185] pumps energy into the dark soliton, which may compensate the 
inhomogeneity-induced emission of radiation, as well as the damping due to the presence 
of the thermal cloud [112j (see Sec. 16.51 below) . 

4-5. Persistence and stability of dark solitons. 

As was highlighted in this Section, there exist many alternative approaches for the 
study of the statics and dynamics of matter- wave dark solitons in the quasi- ID setting. 
Nevertheless, rigorous results concerning the persistence and stability of dark solitons 
in a generalized GP-like model [cf. Eq. (fl8|) ] were obtained only recently [16111186] . 
Particularly, in Ref. |186j . the existence and stability of a black soliton of Eq. ( 1501) were 
studied in the absence of the potential term, while in Ref. |161] a more general model, 
incorporating the potential term, was considered. More specifically, the model used in 
Ref. [16 lj was of the following form, 



where e is a formal small parameter setting the strength of the potential. The results 
obtained in Ref. |161j for bounded and exponentially decaying potentials (as, e.g., ones 
corresponding to red-detuned laser beams — see, e.g., the experiment of Ref. [187] ) can 
be summarized as follows. 

Let us consider that, in the absence of the potential, Eq. f ll06p admits a black 
soliton solution of the form ip(z,t) = q(z — s) exp[— if(n )t + i9] (here, s is the soliton 
center and 9 an arbitrary constant phase), with boundary conditions go ~~ ► =t\/^o as 
z — > ±oo. Then, the above solution persists in the presence of the perturbation induced 
by the potential term in Eq. (11061) provided that the function 



possesses at least one single root, say s - Then, the stability of the dark soliton solution 
depends on the sign of the first derivative of the function in Eq. (I107p . evaluated at s : 



(106) 




(107) 
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an instability occurs, with one imaginary eigenfrequency pair for eM"(sq) < 0, and with 
exactly one complex eigenfrequency quartet for eM"(sq) > 0. In fact, this instability is 
dictated by the translational eigenvalue, which bifurcates from the origin as soon as the 
perturbation is present. For eM"(sq) < 0, the relevant eigenfrequency pair moves along 
the imaginary axis, leading to an instability associated with exponential growth of a 
perturbation along the relevant eigendirection. On the other hand, for eM"(sq) > 0, 
although the eigenfrequency should move along the real axis, it can not do so because 
the latter is filled with continuous spectrum; thus, since the translation mode and the 
eigenmodes of the continuous spectrum have opposite Krein signature, the collision 
of the eigenfrequency of the translational mode with the continuous spectrum results 
in a complex eigenvalue quartet, signalling the presence of an oscillatory instability. 
The relevant eigenfrequencies can be determined by a quadratic characteristic equation 
which, in the case of the cubic GP model (1106!) with f(n) = n and no = 1, takes the 
form |161j . 

A 2 + |m"( So ) =0(e 2 ), (108) 



4 v ' V 2 , 

and the eigenvalues A are related to the eigenfrequencies u through A 2 = — u 2 . For 
sufficiently small e > 0, this equation has only one real root A(e) > for M"(sq) < 
and two complex-conjugate roots, with Re{A(e:)} > for M"{sq) > 0. 

It is interesting to observe that if the characteristic equation (11081) is formally 
applied to the cubic GP model ( 1106D with /(n) = n, no — 1 and V(z) = z 2 , 
one obtains M"(s ) = 2 /_ °° sech 2 (z)dz = 4 and, thus, Eq. (11081) takes the form 
A 2 — (e/2)A + e = 0(e 2 ). Using appropriate rescalings, it can easily be shown that 
the latter characteristic equation can be derived from Eq. fll05f) of Sec. 14.41 Although 
the validity of the radiative boundary conditions for V(z) = z 2 cannot be verified by the 
analysis of Ref. |161j . the above observation leads to the following conjecture |161j : in 
the most general GP model [cf. Eq. fll06|) ] . the two complex-conjugate eigenvalues with 
positive real part for M"(sq) > result from the following Newton's particle equation 
of motion for the soliton center s(t): 

Ho's - e\oM"(s)s = -eM'(s), (109) 

where M(s) is the effective potential implied by Eq. (11071) . while the constants /io and 
A represent, respectively, the soliton's mass and anti-damping — as per the discussion 
of Sec. IH 

The validity of Eq. (I109p . as well as the other theoretical predictions presented in 
this Section, were tested against numerical simulations in Ref. [161] for small decaying 
potentials, and the agreement between the analytical and numerical results was found to 
be very good. Notice that although the above results of Ref. [161] can only be rigorously 
applied to the case of small, bounded and exponentially decaying potentials, the basic 
qualitative features may formally persist for other types of external potentials as well. 
A pertinent example is the work of Ref. [188J, where the the persistence and stability 
of matter-wave black solitons were studied in a condensate characterized by a periodic, 
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piecewise-constant scattering length f : as shown in Ref. |188j . a formal application of 
the predictions of Ref. [161] concerning the persistence and stability of dark solitons in 
this setting, was found to be in very good agreement with relevant numerical findings. 
Nevertheless, an analysis similar to the one presented in Ref. [161] , but for other types 
of potentials (such as confining and periodic ones) is still missing. 



5. Matter-wave dark solitons in higher-dimensional settings 

Quasi-ID matter-wave dark solitons may naturally exist in higher-dimensional settings. 
For example, in the experimentally relevant case of a cigar-shaped trap, the actual 
dimensionality of the BEC density is 3D rather than ID, despite of the fact that the BEC 
can be treated as an effectively ID object using the NLS Eq. ffT8|) with the generalized 
nonlinearities of Eqs. ( !T9|) -(!2T1) — see Sec. 12.41 The density of such a cigar-shaped BEC 
with a dark soliton on top of it is illustrated in the left panel of Fig. [HJ Furthermore, 
quasi-lD dark solitons may also exist in disk-shaped BECs (see Sec. 12. 4p . which are 
described by the (2 + l)-dimensional GP Eq. ( )23l) . The latter, can be expressed in the 
following dimensionless form, 



id t ip 



._ V 2 + V(r) + |# 



0, (HO) 



where V 2 = <9 2 + <9 2 , the density \ip\ 2 , length, time and energy are respectively measured 
in units of 2\/2iraa z , a z , to' 1 and ftw z , while the potential V(r) is given by 

V(r) = iflV, (111) 

with the aspect ratio being Q = u±/u z <C 1. In this case, the soliton of Eq. ( 133|) . 
with the variable z being replaced by x, is an exact analytical solution of Eq. ( II 10p 
for V(z) = 0. This "rectilinear" soliton has the form of a dark "stripe" on top of a 2D 
TF cloud, and the BEC wave function can be expressed (similarly to the ID case) as 
■0 = ipTF( r ) exp(— z/i^^ds^) t)- It is also natural to consider the full (3 + l)-dimensional 
version of Eq. (11 101) . with V 2 = <9 2 + <9 2 + d z , where quasi-lD dark soliton stripes exist 
as well. In this case, the potential and the 3D TF cloud are modified according to 
the relative values of the confining frequencies in the three directions. Examples of the 
densities of a disk-shaped BEC and a spherical BEC carrying a rectilinear dark soliton 
are shown, respectively, in the middle and right panels of Fig. |HJ 

Apart from the quasi-lD dark solitons, purely 2D dark solitons have also been 
predicted to occur in theory (but they have not been observed so far in experiments). 
Such dark soliton solutions of the GP Eq. (jllOp . which have been derived in the 
framework of the small-amplitude approximation (see Sec. 13.41) . may have the form of 
lumps satisfying an effective Kadomtsev-Petviashvili (KP) equation |120] or dromions 
satisfying an effective Davey-Stewartson system |121j ; quasi-lD and 2D dark solitons of 

f BECs with spatially-varying coupling constant g 7 so-called collisionally inhomogeneous condensates 
|189j , have attracted much attention, as they provide a variety of interesting phenomena il 901 4198] . 
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Figure 8. (color online) Examples of the 3D densities of condensates, confined in 
various types of anisotropic harmonic traps and carrying quasi-lD dark solitons. Shown 
are (from left to right) the longitudinal cuts of a cigar-shaped BEC, a disk-shaped BEC, 
and a spherical BEC. 



the dromion type, have also been predicted to occur in disk-shaped mult i- component 
condensates [124]. 

5.1. Snaking instability of rectilinear dark solitons. 

5.1.1. Basic phenomenology and results. An important issue arising in higher- 
dimensional settings is the stability of dark solitons which, for simplicity, will be studied 
at first in the (2 + l)-dimensional geometry (relevant to disk-shaped BECs) and in the 
absence of the potential V(r). The stability of the ID dark soliton stripe (lying, say, 
along the x-direction) in such a 2D setting was first studied in Ref. [199] (see also 
Refs. [20011201] ). In this work, it was shown that the soliton is prone to transverse 
modulational instability, i.e., it is unstable against long- wavelength transverse periodic 
perturbations ~ cos(Qy), where Q is the wave number of the perturbation. In particular, 
the instability band is defined by Q < Q CI , where the critical value of the perturbation 
wave number is given by (for // = !): 



The above expression is a result of a linear stability analysis, which indicates that 
the amplitude of the rectilinear dark soliton will grow exponentially in the transverse 
direction. Nonlinear regimes of this instability were also studied analytically by means 
of asymptotic expansion techniques (see, e.g., Ref. [202] and references therein). This 
instability was extensively studied in the context of nonlinear optics, both theoretically 
[34] and experimentally [203P204] . and was found to be responsible for a possible decay of 
a plane dark soliton into a chain of vortices of opposite topological charges (vortex-anti- 
vortex pairs). Particularly, when the transverse modulational instability sets in, a plane 
black soliton undergoes a transverse "snake" deformation (hence the name "snaking 
instability") [34 ^1202] . causing the nodal plane to decay into vortex pairs. On the other 
hand, unstable gray solitons may not decay into vortices, but rather perform long lived 
oscillations accompanied by emission of radiation in the form of sound waves. 

The basic phenomenology and results described above, persist in the case of other 
higher-dimensional setups as well. For example, in Fig. [9J we show the onset of the 
snaking instability of a rectilinear dark soliton on top of a cigar-shaped BEC (see also 




(112) 
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Figure 9. Contour plots showing the evolution of the density of a cigar-shaped 
BEC confined in the trap V(z,r) = (l/2)(il 2 .r 2 + Q 2 z z 2 ) (parameter values are 
Q, r = 10tt z = 0.1 and \i = 1) and carrying a rectilinear dark soliton. Panel (a) 
shows the initial condition, ip(r,z,0) = yl — V(r, z) tanh(z) and panels (b), (c) and 
(d) are close-up snapshots — at t = 97, t = 101 and £ = 110, respectively — showing 
the onset of the snaking instability and the decay of the soliton into vortex rings. 

Sec. 15.31 below), confined in a trap of the form V(z,r) = (l/2)(Q 2 r 2 + £l 2 z 2 ) (with 
r 2 = x 2 + y 2 ). In such higher-dimensional settings, the unstable soliton collapses into 
more stable vortex structures, namely vortex rings. From the viewpoint of experimental 
observations, the snaking instability and the decay of matter-wave dark solitons into 
vortex rings was first observed in a JILA experiment [48] with a two-component 87 Rb 
BEC (see Sec. I6.ip . In particular, a quasi- ID dark soliton created in one component 
(see Sec. I3.5.3j) evolved in a quasi-spherical trap and, thus, the onset of the snaking 
instability caused the soliton to decay into vortex rings — as predicted in theory [167j . 

5.1.2. Avoiding the snaking instability. As is known from the context of nonlinear 
optics [31], the snaking instability can be avoided by using finite-sized background 
optical beams (see, e.g., relevant experimental results in Ref. [3]). Thus, one should 
expect that the suppression of the snaking instability may also be possible in the case 
of a trapped (disk-shaped) condensate, which also constitutes a background of a finite 
extent. Indeed, in such simple criterion for the suppression of the snaking 
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instability can be found by means of scale competition arguments [21)5] . In particular, if 
the characteristic length scale of the condensate L BEC = 2y/2/Q (i.e., the TF diameter 
for \x = 1) is below the critical length L c = 2ti/Q cx stemming from Eq. (11121) . then the 
snaking instability will not manifest itself. Considering the case of a black (stationary) 
soliton with sin = 0, Eq. ( Ill 2D yields Q cr = 1 and, thus, L c = 2n; in such a case, the 
above scale competition argument, Lbec < L c , leads to the prediction that a use of a 
sufficiently strong trap, such that > f2 c = \/2/n ~ 0.45, can suppress the snaking 
instability. The above criterion was tested against direct numerical simulations [205J, 
and it was found that the critical value of the trap strength is less than the theoretically 
predicted, namely Q c ~ 0.31. This discrepancy can be understood by the fact that for 
small BECs (i.e., for tight traps) the presence of the dark soliton significantly modifies 
the maximum density which is less than /i by a "rescaling" factor /, found to be / ~ 0.5. 

On the other hand, it was recently predicted 2()(i that stable 3D stationary 
dark solitons may exist in dipolar condensates (for this type of BECs see, e.g., the 
recent review |207] and references therein). In particular, the special feature of dipolar 
condensates, namely the dipole-dipole interaction, together with the use of a sufficiently 
deep optical lattice in the soliton's nodal plane, allows for the existence of dark solitons 
of arbitrarily large transversal sizes, which are not prone to the snaking instability. 
In this case, the underlying reason for the suppression of the snaking instability is 
that the dipole-dipole interaction is long-range (it decays like r -3 , where r is the inter- 
particle distance), which means that the respective GP equation incorporates a nonlocal 
nonlinear term. Generally, such a nonlocal response may arrest collapse and stabilize 
solitons in higher-dimensions, as was shown in the context of optics (see, e.g., Ref. [208J, 
as well as Ref. [209J for a relevant recent work on dark solitons). 

We also note that a more "exotic" dark soliton configuration in the 2D setting, 
which is not subject to the snaking instability, was presented in Ref. |210j (see also 
Ref. [21 lj ) . This configuration, which refers to a two-component BEC (see Sec. 16. ip . 
consists of a "cross" formed by the intersection of two rectilinear domain walls, with 
the wave functions of the same species filling each pair of opposite quadrants having a 
7r phase difference. This way, a dark soliton configuration is formed, which was found 
to be stable for long times — and even in the presence of rotation of the trap — in a 
large parametric region. 

5.2. Matter-wave dark solitons of radial symmetry 

5.2.1. Ring dark solitons (RDS) and spherical shell solitons (SSS). A special class 
of dark solitons in higher- dimensional settings consists of dark solitons exhibiting a 
radial symmetry, that can be realized by wrapping the nodal plane around on itself. 
Such structures were first introduced in the context of nonlinear optics |212j . with the 
motivation being that these dark solitons may not be prone to the snaking instability: 
indeed, if a dark stripe is bent so as to form a dark ring of length L < 2n/Q CT (in the 
(2 + l)-dimensional geometry), then the snaking instability will be suppressed. Such 
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Figure 10. (color online) Examples of the 3D densities of condensates, confined in a 
disk-shaped (left panel) and a spherical (right panel) traps, and carrying a ring dark 
soliton and a spherical shell soliton, respectively. 



ring dark solitons (RDSs) were studied in nonlinear optics both theoretically [21311214] 
and experimentally [21514217] . while later were also predicted to occur in BECs [218J. In 
this context, and in the 2D setting (i.e., in a disk-shaped BEC), the RDS has the form 
of an annular "trough". On the other hand, in a 3D setting (i.e., in a spherical BEC), 
the radially symmetric dark soliton is called spherical dark soliton [158] . or spherical 
shell soliton (SSS) — according to the nomenclature of Ref. [219] — and has the form of 
a nodal spherical "shell" . Examples of the 3D densities of a disk-shaped and a spherical 
BEC, carrying a RDS and a SSS, are shown in Fig. [TUJ 

As was originally proposed in Ref. [218] . RDS may be generated in BECs by 
means of phase-engineering techniques (i.e., by a proper phase-imprinting method — see 
Sec. l3.5.1l below). similar to the ones used for the generation of optical RDS [215H217] . 
Another technique that has been proposed for the generation of RDS in BEC is the 
matter- wave interference method (see Sec. 16. 2j) : if the condensate is initially trapped in 
a narrow cylindrical box-like potential, and then is allowed to coherently expand in the 
presence of a wider cylindrical impenetrable hard-wall potential, it is reflected from the 
boundary, and the self-interference pattern has the form of a sequence of non-stationary 
concentric RDS [2201- 1222] (see also relevant work in Ref. [223]). 

5.2.2. Dynamics and stability of RDS and SSS. From a mathematical standpoint, 
matter-wave dark solitons of radial symmetry can be considered as quasi-lD objects 
and, accordingly, be analyzed by means of a quasi-lD GP equation. In particular, 
either RDS or SSS can be described by Eq. (jllOp . with the Laplacian being in the form, 

V 2 = ^ + ^^«9 r , (113) 

r 

with D — 1, 2, 3. In this setup, the simplest case of D = 1 reduces Eq. f II 1 j) to the ID 
GP Eq. (T2"4"|) describing a quasi-lD BEC (here, r = z). The higher- dimensional setups 
correspond to the cases of D = 2 and D = 3: in the former case, Eq. (IllOp describes a 
disk-shaped BEC in the (x, y) plane (with r given by r = a/x 2 + y 2 ), while in the latter 
case Eq. (11101) describes a spherical BEC (with r given by r = a/x 2 + y 2 + z 2 ). 
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In such a quasi-lD setup, the Hamiltonian or the Lagrangian perturbation theory 
for dark solitons (see Sec. 14. 2[) may also be applied for the study of the dynamics 
of RDS and SSS. In particular, Eq. ( 11 101) . with the Laplacian of Eq. (I113p . can be 
treated as a perturbed ID NLS equation [similar to Eq. (!7H|) ] provided that both the 
potential term and the term ~ r _1 can be considered as small perturbations. The 
latter assumption is physically relevant for radially symmetric solitons of large radius 
r . Then, approximating the functional form of the RDS or SSS as [cf. Eq. (186|) ]. 

ip s (r, t) = cos 4>{t) tanh ( + i sin 4>{t), ( = cos 4>(t) [r — r (t)] , (H4) 

where r (t) is the soliton radius, it can be found |218] (see also Ref. |158j ) that r is 
governed by the following Newtonian equation of motion: 

d 2 r l dV e s . . 

Here, the effective potential is given by 

V eS (r ) = V(r )-lnr 2 (D ~ 1)/3 , (116) 

and V(r ) = (l/2)Q 2 rQ is the trapping potential evaluated at the soliton radius r . In 
the ID limit of D = 1, the last term in the right-hand side of Eq. ( 11161) vanishes and 
Eq. (I115p is reduced to Eq. (1901) . In the higher-dimensional cases of D = 2 or D = 3, 
the equation of the soliton motion ( 1115P is clearly nonlinear (even for nearly black RDS 
or SSS) due to the presence of the repulsive curvature-induced logarithmic potential. 

Equation (11151) predicts the existence of both oscillating (grey) and stationary 
(black) RDS or SSS: the former, perform oscillations on top of the TF cloud, changing 
their radii between a minimum and a maximum value, while the latter correspond to 
the minimum of the effective potential Eq. ( 11161) (such stationary states do not exist in 
the case of a uniform ground state — as in the context of nonlinear optics [212] ). As 
was shown in Ref. [218J (see also Ref. |158j ). such a particle-like approach can describe 
quite effectively the above generic scenarios and the RDS dynamics up to a certain time 
(before the development of instabilities — see below). Furthermore, in Refs. [22411225] it 
was shown that the dynamics of small-amplitude RDS, as well as the collisions between 
them, can be described in the framework of an effective cylindrical KdV equation |114] 
(see also Refs. [212T l2l4j for similar findings in optics). On the other hand, numerical 
simulations of Ref. |218] revealed that RDS are generally unstable, as they either decay 
to radiation (the small-amplitude ones) or are subject to the snaking instability (the 
moderate- and large-amplitude ones). Interestingly, as shown in Fig. [HJ the snaking 
instability of the RDS results in the formation of vortex-anti- vortex pairs in multiples of 
four, which are initially set along a ring, forming a so-called vortex necklace. Eventually, 
this pattern relaxes to a set of four pairs located on a ring, which oscillates in the radial 
direction between the same limits which confined the oscillations of the original RDS; 
simultaneously, the pairs perform oscillatory motion along the ring [218J. 

Matter-wave dark solitons of radial symmetry were also analyzed by means of other 
approaches. For example, in Ref. |219] (see also Ch. 7 in Ref. jl3]) RDS and SSS were 
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Figure 11. (color online) Contour plots of the density of a disk-shaped condensate 
carrying a stationary RDS, which develops the snaking instability. The initial condition 
used for the integration of the GP Eq. (|110j) is ->p(r,0) = \i - (l/2)Q 2 r 2 tanh(r - 
r (0)), with n=l,Sl = 0.035 and initial soliton radius r (0) = (V2D,)- 1 = 20.2 (for 
a discussion concerning the value of r (0) for stationary RDS see Ref. [158 ). The top 
left (top right) panels show, respectively, the initial condition (t = 0) and a snapshot at 
t = 40, while the bottom left (bottom right) right panels show, respectively, the onset 
of the snaking instability (t = 80) and the formation of the vortex necklace (t — 100). 



considered as nonlinear Bessel functions, namely solutions of the equation, 

q" + -q' ~^q + 2/ig - 2g 3 = 0, (117) 

resulting from Eq. ( Ill OR (with V(r) = 0) when the ansatz ip = q(r) exp(— i/it + iScp) 
is introduced (in the latter expression, S is the topological charge of a central vortex). 
In this setting (i.e., in the absence of the trap), the solutions of Eq. ( 11171) include, 
apart from the ground state, singly- and multiply-charged vortices, as well as infinitely 
many RDS; the nodes of the latter correspond to the nodes of the nonlinear Bessel 
function governed by Eq. (11171) . On the other hand, if an external harmonic potential 
is incorporated in Eq. ( 11171) . then it is possible to find infinite branches of nonlinear 
bound states, with each branch stemming from the respective mode of the underlying 
linear problem (the radially-symmetric 2D quantum harmonic oscillator) [226] . 

A stability analysis performed in Refs. [219(1226] also revealed that the radially- 
symmetric dark solitons are typically unstable but, in agreement to the findings of 
Ref. |218j . their life-times may be considerably long. Additionally, as shown in a recent 
work [227] . the life-time of RDS may be extended employing the so-called [228] Feshbach 
Resonance Management (FRM) technique, which is based on the use of external fields 
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to periodically modulate in time the s-wave scattering length [229- 1233] . In any case, 
the theoretical investigations indicate that RDS and SSS have a good chance to be 
experimentally observed. In fact, structures similar to stationary SSS have already 
been observed as transients in the Harvard experiment of Ref. [65J (see also Ref. [234J 
were SSS are predicted to occur as a result of collisions of vortex rings). 

5.3. Stability of dark solitons in cigar-shaped condensates. 

The transverse (in) stability of dark solitons confined in a purely 3D setting, namely in 
a cylindrical trap of the form V(z,r) = (l/2)m[u 2 z 2 + uj_(x 2 + y 2 )}, was first studied 
in Ref. |145] . In that work, it was shown that dynamical stability of a black soliton 
(stationary kink), say ipo, with a nodal plane perpendicular to the axis of the cylindrical 
trap (see left panel of Fig. E}, requires a strong radial confinement [as in the (2 + 1)- 
dimensional case]. Particularly, it was shown that the instability can be suppressed if the 
transverse (radial) condensate component is not in the TF regime, which is guaranteed 
as long as ftw± > [l (where [i is the 3D chemical potential). This criterion can physically 
be understood as follows: if the condensate is confined in a highly anisotropic (cigar- 
shaped) trap, then the energy of the lowest possible radial excitation, Jtw±, must exceed 
the kink-related kinetic energy, K = —(l/2)ip d 2 ipo ~ so that the latter can not 
be transferred to the BEC's unstable transverse modes by the inter-atomic interaction. 
Moreover, a systematic study in Ref. [145] revealed a criterion of dynamical stability for 
the black soliton in terms of the ratio u±/u z , namely, 



The critical value 7 C was calculated for various values of u±/u z and it was found that a 
characteristic value of 7 C , pertinent to the limiting case of u± ^> u z , is 7 C ~ 2.4. 

In a more quantitative picture, a detailed study of the BdG equations in Ref. [145] 
(see also relevant work in Refs. |167yi68j ) revealed the emergence of complex eigenvalues 
in the excitation spectrum and their connection to oscillatory dynamical instabilities, 
including the snaking instability. Additionally, in Ref. |235] it was found that the 
emergence of complex eigenvalues is directly connected to bifurcations of the rectilinear 
black solitons to other stationary states that may exist in BECs confined in cylindrical 
traps. In particular, an investigation of the dependence of the excitation energy on the 
dimensionality parameter d [cf. Eq. (1161) ] led to the following results: for sufficiently 
low excitation energies, the black soliton may bifurcate to a solitonic vortex or to 
an axisymmetric vortex ring (see also Ch. 7 in [43] and references therein), with the 
corresponding bifurcation points occurring at a low and a higher value of d. It was also 
found that the emergence of the first (second) complex eigenvalue in the Bogoliubov 
spectrum coincides with the above mentioned bifurcation points. Therefore, according 
to the above results, it can be concluded that the emergence of the complex eigenvalues 
in the excitation spectrum (a) denote the onset of dynamical instabilities of black solitons 
and (b) indicate the excitation of lower energy topological states. Since these states are 
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energetically preferable, the onset of the dynamical instability will result in decay of the 
"high-energy" black soliton to these "low-energy" states carrying vorticity. 

On the other hand, the stability of moving (gray) solitons was analyzed in Ref. |91j . 
According to this work, and following the arguments of Ref. |145j . a criterion for not 
being in the radial TF regime (which is required for dynamical stability of the gray 
soliton) is w < R ~ £, where R is the radial size of the BEC and w is the soliton width 
(£ is the healing length). Recalling that the soliton width w and velocity v depend on 
the soliton phase angle as w ~ 1/ cos0 and v ~ sin0, it is clear that as the soliton is 
moving towards the boundaries of the BEC, its width (velocity) is increased (decreased). 
Thus, the instability border R ~ w for the gray soliton is reached for larger values of 
the parameter 7 [see Eq. f U18p ] than the ones pertaining to the black soliton. In the 
case of strongly anisotropic traps, the critical value of the chemical potential required 
for the dynamical stability of the gray soliton is proportional to the soliton velocity. In 
other words, the stability domain of gray solitons is wider than the one of black solitons. 
In fact, the shallower the soliton it is, it becomes more stable, similarly to the case of 
homogeneous systems: see Eq. (I112p . which indicates that the instability band vanishes 
for shallow solitons with cos0 — > 0. 

Numerical simulations of Ref. [91] have also revealed that for 7 > 10, a phase- 
imprinted dark soliton — with a 7r-phase jump — always decays in a cigar-shaped BEC 
(on a time scale of order of u/[ ), while for 7 < 5 it transforms into a dark soliton 
characterized by a flat notch region and r-independent velocity. Here, it is relevant to 
mention that in the recent Technion experiment [70] , an interesting nonlinear excitation 
that evolves periodically between a dark soliton and a vortex-ring was observed in a 
87 Rb BEC, for 7 = 4.95 (see also relevant theoretical work in Refs. [23611237] ). 

5.4- Matter-wave dark solitons in the dimensionality crossover from 3D to ID 

5.4-1- The single- soliton state. As discussed in Sec. l2.4[ if the dimensionality parameter 
[cf. Eq. ( fl6l) ] takes values d ~ 1, then a cigar-shaped condensate is in the so-called 
dimensionality crossover regime from 3D to ID. However, in such an experimentally 
relevant regime f , exact analytical dark soliton solutions (of arbitrary amplitudes) of the 
pertinent effectively ID mean-field models (see Sec. 12.41) are not available. As a result, 
the analytical techniques presented in Sec. 14. 21 cannot be applied for the study of matter- 
wave dark soliton dynamics in this regime. Nevertheless, the results of Sees. 14.21 and 15.31 
indicate that the evolution of dark solitons should be similar to the one pertaining to 
the TF-1D regime, while the solitons would not be prone to the snaking instability, as 
in the case of the purely higher-dimensional setups. 

The statics and dynamics of matter-wave dark solitons in the crossover regime 
between ID and 3D were studied in Ref. [95]. Particularly, numerical simulations of 
the 3D GP equation revealed that matter-wave dark solitons are indeed dynamically 
stable, and perform harmonic oscillations in the harmonic trap. Importantly, in the 

f Note that the recent Heidelberg experiments [69l[TT] were conducted in this regime. 
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case of small- amplitude oscillations, the oscillation frequency oo osc resulting from the 3D 
GP equation was found to be equal to the eigenfrequency oja of the anomalous mode 
of the effectively- ID NPSE model [cf. Eqs. ( fl8l) and ( fl9l) ]; note that the latter, can be 
expressed in the dimensionless form: 



id t ip 



A (H9) 



2 2 2 2(1 + |V| 2 ) 1/2 . 

where units are the same to the ones used for Eq. f[2~4l) . The above finding leads, in turn, 
to the following conclusion: an equation of motion for the center zq of a dark soliton in 
a condensate in the dimensionality crossover regime can be expressed as follows: 

d 2 z dV cS 1 2 2 _ 

IF = = 2 W ° scZ °' Uosc = Ua ' ^ ' 

As shown in Ref. [OS], the soliton oscillation frequency is a decreasing function of the 
dimensionality parameter d, taking values ranging from u osc = Q (corresponding to the 
non-interacting limit of d — > 0), to tu osc = Q/y/2 [corresponding to the TF-1D regime of 
d <C 1 — cf. Eq. (19~T]) ]. with Q being the normalized strength of the harmonic trap. In 
any case, the oscillation frequency is up-shifted from its TF-1D value and, as a result, 
the effective trapping potential felt by the dark soliton during its motion [see Eq. ( MI]) ] 
will effectively become steeper. In that regard, it is relevant to mention that substantial 
shifts of the oscillation frequency (which may be of order of 10%) have been predicted 
in Ref. [95] and later confirmed in the Heidelberg experiment of Ref. [69J (see discussion 
in Sec. I5.4.3p . It should also be noticed that the soliton oscillation frequency is also 
a decreasing function of the soliton amplitude (or, in other words, of the oscillation 
amplitude), contrary to the result corresponding to the TF-1D regime f [69|[TTj . 

Results similar to the ones obtained in the framework of the NPSE model [95] can 
also be obtained in the framework of the generalized NLS Eq. ([18]) with the nonlinearity 
function of Eq. (120]) (see, e.g., the analysis of Ref. [71]). Furthermore, we should mention 
that the oscillations of dark solitons were also analyzed in the framework of a GP model 
with generalized nonlinearities, and specific results in the physically relevant case of a 
cubicquintic nonlinearity (modeling two- and three-body interactions — see Sec. 12.41) 
were presented |162j . The same model was also studied in Ref. [238] , were various 
stationary states, including dark solitons, were found and analyzed in detail. 



5.4-2. The multiple- soliton state. Apart from the single-dark soliton state, the case 
of a multiple-dark soliton state can also be analyzed in the dimensionality crossover 
regime using, as a guideline, the methodology exposed in Sees. 13.61 14.21 and 14.31 In 
particular, in the case of well-separated and symmetrically interacting dark solitons in 
a harmonic trap, one may follow the analysis of Refs. J69J[7T] and analyze this problem 
by adopting a simple physical picture: each soliton in the multi-soliton state follows 
the evolution of the single-soliton state, i.e., it oscillates in the trap with an oscillation 

f Recall that the results of Sec. 14.21 indicate that the soliton oscillation frequency does not depend on 
the soliton amplitude in the TF-1D regime. 
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Figure 12. (Color online) The lowest characteristic eigenfrequencies of the Bogoliubov 
excitation spectrum for a condensate, in the dimensionality crossover regime from 3D 
to ID, carrying two dark solitons. The model used is the NPSE (I119p . with parameters 
£1 = 0.05, and \l = 1.86; the dimensionality parameter is d — N£la/a± = 0.82. The 
lowest characteristic eigenfrequencies of the Bogoliubov excitation spectrum: shown 
are the one located at the origin, the one at O = 0.05, the one at \^3Q = 0.087 
(corresponding to the Goldstone mode, the Kohn mode and the quadrupole mode, 
respectively), as well as the two anomalous modes, one with u>i = 0.756J1 = 0.0378 
and one with uj 2 = 2.094^ = 0.1047. 



frequency u osc equal to the eigenfrequency u)\ of the first anomalous mode (determined 
by a BdG analysis of the pertinent ID mean- field models of Sec. I2.4p . and interacts with 
the neighboring solitons via the effective repulsive potential of Eq. ( 1761) . 

In order to further elaborate on the above, let us consider — as an example - 
the two-dark soliton state of the NPSE model of Eq. f II 1 9 j) . This state (which, in 
the linear limit, corresponds to the second excited state of the quantum harmonic 
oscillator) can be obtained as a nonlinear stationary state of the system by means 
of, e.g., a Newton- Raphson method. The pertinent configuration, has the form of two 
overlapping dark solitons, placed at Zi = ±2.185 (i = 1, 2) with a fixed relative distance 
5zq = 4.37. The corresponding Bogoliubov excitation spectrum, namely the spectral 
plane (cu r , Uj) of the eigenfrequencies oj = u r + iui, is shown in Fig. [12j as it can clearly 
be observed, among the lowest eigenfrequencies (such as the ones at uj = 0, Ud = Q and 
co q ~ y/3Vl, corresponding to the Goldstone, dipole and quadrupole mode, respectively), 
there exist two anomalous modes with eigenfrequencies oj\ = 0.75QQ = 0.0378 and 
u 2 = 2.094ft = 0.1047. 

According to the analysis of Sec. I4.3[ small displacements of the dark solitons from 
their equilibrium points lead to the in-phase and out-of-phase oscillatory motion of the 
dark soliton pair (see Fig. [13]), with the respective oscillation frequencies being equal 
to the eigenfrequencies U\ and u 2 of the two anomalous modes. Importantly, the value 
of the eigenfrequency of the first anomalous mode, u\ = 0.0378, is quite close to the 
oscillation frequency of a single dark soliton, o; osc = 0.0375, in the same setup (i.e., with 
the same parameter values), with the percentage difference being w 1.6%. This generic 
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example suggests that, generally, the dynamics of the two dark soliton state can be 
described by an effective Lagrangian for the two solitons, namely L e g = T — V; here, 
T and V are the kinetic and potential energies, respectively, depending on the soliton 
centers, Zi (i = 1,2), and soliton velocities, = d,Zi/dt, as follows, 

2 2 

T = E^) 2 > V = J2^u 2 osc z^ + V hlt (z 2 -z 1 ), (121) 

i=l i=l 

where V int (z 2 — Z\) (with z 2 — Z\ = 2z ) is the repulsive potential of Eq. fl76|) . Then, 
the evolution of the soliton centers can readily be determined by the Euler-Lagrange 
equations d^d^L^) / dt — d Zi L cS = 0. The latter, may be simplified upon using the 
approximate form of the repulsive potential [cf. Eq. f!77p ]. thus leading to the following 
equations of motion: 

zi= - ul^zx - 8n 3 /2 exp[-2y/n^(z 2 - z{)], (122) 
z 2 = - u 2 osc z 2 + 8n 3 /2 exp[-2y/n^(z 2 - z\)], (123) 

where we have assumed well-separated, almost black solitons [i.e., in Eq. ( 1771) we have 
set B pd 1]. Apparently, the above analysis can readily be generalized for multiple 
solitons [7_I], with each one interacting with its neighbors. Importantly, if u osc in the 
system of Eqs. (I122[) ( TT23T) was considered to be unknown, then it would be possible to 
be directly obtained in the form of the characteristic frequencies of the normal modes 
of this system (see details in Ref. [H]); these characteristic frequencies coincide to the 
ones determined via the BdG analysis. This results justifies a posteriori the considered 
decomposition of the principal physical mechanisms (oscillations and interactions of 
solitons) characterizing the system. 

We should also note that apart from the case of small-amplitude oscillations of 
two well separated, almost black solitons, the more general case of the dynamics of 
n-interacting dark solitons (which may also perform large amplitude oscillations — see 
Sec. 15.4.31 below) is possible using the full set of the above mentioned Euler-Lagrange 
equations; the latter, lead to the following n-coupled equations of motion [71] : 

>A / d 2 V d 2 V \ dV , x 

k=l 

where V = Yn=i v i is tli e potential energy, V { = n oB 2 j/ {2 sinh 2 [ y /n^B ij (z i -z j )}} is 
the interaction potential felt by the i-th soliton due to the presence of the other solitons, 
while Zij = (1/2) (zi — Zj) and Bij = (l/2)(Bj + Bj) denote, respectively, the relative 
coordinate and the average depth for solitons i and j. 

5.4-3. Large amplitude oscillations and experimental observations. Let us now return 
to the above example of the two-dark soliton state of Eq. (11191) and consider again 
the out-of-phase oscillation (we will use the parameter values of Fig. fl~2l) . In this 
case, if the initial soliton separation is significantly larger than 5z$ = 4.37 — or, in 
other words, if the displacements of solitons around their equilibrium positions are not 
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Figure 13. (Color online) Top panel: Spatio-temporal contour plot of the density of 
a cigar-shaped BEC confined in a trap of strength = 0.05. The condensate is in 
the dimensionality crossover regime from 3D to ID, and the model used is the NPSE 
Eq. (|119[) . The dark solitons, initially placed at z = ±4 (i.e., Sz = 8 > Szo — 4.37), 
oscillate out-of-phase with a frequency lu osc = 1.740 < ui2 = 2.09417. Bottom panel: 
The oscillation frequency of the dark solitons, w osc , as a function of their initial relative 
distance Sz. For Sz — » Szq (corresponding to the stationary state), we obtain w osc — > LU2, 
while for Sz ^> Szq, we obtain w OS c — ► 2wi. 

small — then u osc differs from (in fact, it is quite smaller than) u^: considering, e.g., 
that 5z = 8 (corresponding to initial locations of the soliton centers z = ±4 — see 
Fig. [13]) , the out-of-phase oscillation of the two solitons is characterized by a frequency 
u osc = 1.74Q < uj2 = 2.09AQ. A qualitative explanation for this difference is the 
following: As discussed above, the evolution of two initially overlapping dark solitons 
can be effectively described by the equations of motion (I122p - (1123p . in the presence of 
the repulsive potential which depends exponentially on their relative distance. If the 
relative distance between the two solitons is not significantly different than the one 
pertaining to the corresponding stationary state, i.e., 5z w 5z , the effective repulsive 
force is strong and, as a result, their motion is strongly affected by their coupling. On 
the other hand, if their initial separation becomes larger (as, e.g., in the case of the 
example under consideration, with 5z = 8), the repulsive force becomes exponentially 
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small and, as a result, the motion of each individual soliton is not significantly affected 
by the presence of the other. 

In the bottom panel of Fig. [13] we show the oscillation frequency of the dark soliton 
- obtained by the direct integration of the NPSE model of Eq. (11191) — as a function 
of the initial relative distance between the two solitons. It is clear that the oscillation 
frequency, which takes values in the interval 5z > 5zq (the value corresponding to 
the stationary state), exhibits two different asymptotic regimes: when the initial soliton 
separation is small enough, 5z — > 5zq (i.e., for strong coupling between the two solitons), 
the oscillation frequency tends to the eigenfrequency Co> 2 of the largest anomalous mode; 
on the other hand, when the initial soliton separation is large enough, 5z ^> 5zq (i.e., 
when the solitons are actually decoupled), the oscillation frequency tends to %jJ\\ the 
latter value can be explained by the fact that the period of oscillation of each individual 
soliton is the half of the one that would correspond to a single soliton oscillation. 

As concerns relevant experiments, an oscillating and interacting dark soliton pair 
in a 87 Rb BEC, in the dimensionality crossover regime between ID and 3D, was 
experimentally observed in the Heidelberg experiment of Ref. [69]. In this experiment, 
large amplitude oscillations were induced by the method of matter-wave interference (see 
Sec. !6.2l below). The dependence of the oscillation frequency on the distance between the 
two solitons (see the left panel of Fig. ITS]) was found and compared with experimental 
data: the agreement between the theoretical predictions (based on a study of the 
NPSE model) and the experimentally observed oscillation frequencies was excellent. 
In accordance to the analysis of this Section, considerable upshifts — up to 16% — of 
the soliton oscillation frequency from the value of Q/y2 were observed in the study of 
Ref. [69], and they were quantitatively attributed to the dimensionality of the system 
and the soliton interactions. 

6. Matter-wave dark solitons in various settings and parameter regimes 

6.1. Matter-wave dark solitons in multi- component condensates. 

Multi-component ultracold atomic gases and BECs may be composed by two or more 
atomic gases, which may have the form of mixtures of: (a) two different spin states 
of the same atom species (so-called pseudo-spinor condensates) [239-241J; (b) different 
Zeeman sub- levels of the same hyperfine level (so-called spinor condensates) |242] ; (c) 
different atom species (so called heteronuclear mixtures) |243] : (d) degenerate boson- 
fermion clouds |244] : (e) purely degenerate fermion clouds |245] (see also [4TH43] for 
reviews and references therein). Such multi-component systems support various types 
of matter-wave soliton complexes, with the type of soliton in one species being the same 
or different to the type of soliton in the other species. Here, of particular interest are 
the so-called vector solitons with the one component being a dark soliton, which have 
mainly been studied in the context of two-component and spinor condensates. 
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6.1.1. Dark solitons in two- component condensates. Generally speaking, a mixture of 
M purely bosonic components can be described in the framework of mean-field theory 
by a system of Af coupled GP equations, which can be expressed in the following 
dimensionless form: 

dib 1 ^ 

i-QjT- = ~2 V2 ^n + V n (r)lp n + ^ [9n,k\*Pk\ 2 lpn ~ ^n,k^k + ^n,k^n] ■ (125) 

k=l 

Here, ip n is the wave function of the n-th component (n = l,...,Af), V n {r) is 
the trapping potential confining the n-th component, A n ^ is the chemical potential 
difference between components n and k, the nonlinearity coefficients g n ^ = 9k,n 
characterize inter-atomic collisions, while the linear coupling coefficients = Kk,n 
account for spin state inter-conversion, usually induced by a spin-flipping resonant 
electromagnetic wave (see, e.g., Ref. |246j ). In some works (see, e.g., Refs. [247H249] ) 
fermionic mixtures are also described in the framework of the mean-field theory, with 
the self-interacting nonlinear terms being replaced by g n ,n|V'n| 4/ ' 3 Vv Notice that in the 
GP Eqs. (I125j) both the energy E and the total number of atoms, N = ^2 k=1 Nf. = 
Y^,k=i I lV'fc| 2 ^ r ) are conserved; furthermore, in the absence of linear inter-conversions 
(Kn,k — 0), the number of atoms of each component N k is separately conserved. 

Let us consider the case of two bosonic species (Af = 2), and assume that the 
system is homogeneous (V n = 0). If, additionally, there is no spin-state inter-conversion 
(«n,fc = 0) and chemical potential difference = 0), then the binary mixture is 

immiscible provided that the following immiscibility condition holds [250J, 

A = (g\292\ - 9xx922) /9n > °> ( 126 ) 

where the, so-called, miscibility parameter A takes in practice the values A?s9x 10 -4 
or A rj 0.036 for a mixture of two spin states of a 87 Rb BEC [2391 EMU] (see also 
Ref. [262]) or a 23 Na BEC [21], respectively. The condition (1X261) indicates that if the 
mutual repulsion between species is stronger than the repulsion between atoms of the 
same species then the two species do not mix. In such a case, the two species tend 
to separate by filling two different spatial regions, thus forming a "ball" and "shell" 
configuration (see, e.g., Ref. [240] for relevant experimental results). This way, the 
ground state of the system — i.e., the state minimizing the energy — may take the 
form of domain-wall solutions of the GP Eqs. f 1 1 2 5 j) [251H256] . In accordance to the 
experimental observations, these solutions represent configurations of the following form: 
in the Thomas- Fermi limit (where kinetic energy is negligible), one species occupies the 
region around the trap center, and it is separated by two domain- walls from side domains 
occupied by the other species; on the other hand, kinetic energy favors a configuration 
where a single domain-wall at the trap center separates two domains occupied by the 
different species [255] . The dynamics of phase-separation of two-component BECs has 
been studied in various works both theoretically (see, e.g., Refs. [25711259] and also 
Refs. [26011261] for proposed applications) and experimentally [262H264] . Importantly, 
magnetic-field Feshbach resonances can be used to controllably change the the inter- 
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species [263] or the intra-species |264] scattering length, and thus controllably change 
the (im)miscibility between the two species [264J. 

Apart from domain- walls, a trapped two-component quasi- ID BEC supports vector 
solitons, with the one component being a dark soliton; in such a case, typically, the 
other component may be a dark soliton |122U124l 12581 12651 - 1271] or a bright soliton 
[131[ 12581 1269[ [272-274] . Additionally, apart from such dark-dark and dark-bright 
solitons, dark-antidark f solitons have also been predicted to exist in BECs, either 
in a stationary form [267] or as a dynamical entity resulting from instabilities [258, 270J . 
Below we will focus on the most generic vector matter-wave solitons, namely the dark- 
dark and dark-bright ones (note that the latter have also been observed experimentally 
[67J), presenting results corresponding to the simplest possible setup, which solely 
includes the traditional time-independent harmonic trap. Notice that in the absence 
of the trap, vector solitons of the above mentioned types have been extensively studied 
in the context of nonlinear optics: there, multi-component solitons occur when fields 
of one frequency, or one polarization, become coupled to fields of other frequencies, or 
other polarizations (see, e.g., the review [34J and references therein). Mathematically 
speaking, the existence (and stability) of multi-component optical solitons (and also 
matter- wave solitons in the miscible case) can be understood by the fact that the relevant 
coupled NLS equations rely on the so-called Manakov system |275] : the latter, has the 
form of a vector NLS equation, namely, 



where u(z,t) = (ui(z,t),U2(z,t), ■ ■ ■ ,u n (z,t)) is a n-component vector. This system is 
known to be completely integrable [27611278] (in fact, it can be integrated by extending 
the 1ST method that has been used to integrate the scalar NLS equation [H1IT2]) and 
admits such vector iV-soliton solutions [279H281] . 

As shown in Refs. [265,266], the dynamics and interaction of dark-dark solitons 
in a two-component quasi-lD BEC can be studied by means of a variational approach; 
in the case of equal chemical potentials, the latter is based on the use of the following 
ansatz for the single- component soliton wave-functions, 



where 2z (t) denotes the relative distance between the two solitons, and the =F signs 
correspond, respectively, to a kink-anti-kink state (where the solitons' phase fronts are 
facing each other) and a kink-kink state (where the solitons' phase fronts are in the same 
direction). Both the miscible and immiscible cases where studied in Refs. |265p fi6j and 
the main results of the analysis can be summarized as follows. In the miscible case 
(#11 = 922 = 912)1 and for the kink-anti-kink state, the trajectories in the (zq, zq) phase- 
plane are either periodic surrounding the center (0, 0) [indicating the formation of a 

f An anti-dark soliton is actually a dark soliton with reverse-sign amplitude, i.e., it has the form of a 
hump (instead of a dip) on top of the background density (see, e.g., Refs. [117H118] ). 



id t u 




u ± |u| 2 u, 



(127) 



i[>i(z,t) = 5tanh [B {z - z (t))}+ iA, 
ipz(z,t) = 5tanh [B (z + z (t))] =F iA, 



(128) 
(129) 



59 



bound state ("soliton molecule")], or free [indicating acceleration (deceleration) of the 
approaching (outgoing) solitons]; contrary, in the case of the kink- kink state, where 
solitons move in the same direction, the solitons form a bound state which can never 
be broken. On the other hand, in the immiscible case (i.e., when domain- walls are 
present), it was shown that if a dark soliton exceeds a critical velocity then it can be 
transferred from one component to the other at the domain-wall; on the other hand, for 
lower velocities, multiple reflections within the domain were observed. In such a case, 
the soliton is accelerated after each reflection and eventually escapes from the domain. 

As mentioned above, dark-bright matter-wave solitons in a quasi-lD binary BEC 
are also possible. Particularly, in the miscible case (with all nonlinearity coefficients 
g n> k being normalized to unity), the wave functions ipd( z , t) and ipb(z, t) of the dark and 
bright soliton components may be expressed in the following form |131j . 

ipd(z, t) = A/zIcos^tanh {k [z — -Zo(t)]} + iy/J^sincp, (130) 
ipb{z,t) = \J ^N h K sech {k [z - z (t)}} exp(z0 b ) . (131) 

Here, = fi and /x b = fi + A are the chemical potentials of the dark and 
bright components, <fi is the dark soliton's phase angle, zq denotes the solitons' 
center, N\> = \iph\ 2 dz is the normalized number of atoms of the bright soliton, 
k = a/ Jl cos 2 (j) + (iV b /4) 2 — iV b /4 is the inverse width of the bright soliton, and 
#b = (fttan0)x+ [ft 2 (l — tan 2 (f>)/2 — A]t is the bright soliton's phase. According to the 
analysis of Ref . |131] , if the external trapping potentials Vd and for the dark and bright 
solitons are slowly varying on the soliton scale K~ l , then the dynamics of the dark-bright 
soliton can be described by the effective particle approach of Sec. 14.21 Particularly, 
assuming that the solitons are sufficiently slow, a multiple-time-scale boundary-layer 
theory — similar to the one used in Ref. |144] - leads to the following equation of 
motion for the soliton center, 

^ = -i^o) - »o)-« (132) 
dt 2 2 dl 0) 8 v //i+(iV b /4) 2 -K 1 (z )' 1 J 

where V^ h {z ) = dV A , h /dz Q . In the limit iVb -> 0, Eq. ffT32l) is reduced to Eq. (1901) 
(recall that the latter predicts dark soliton oscillations with a frequency VL/y/2), while 
the motion of the vector soliton becomes more sensitive to the presence of the bright 
component as is increased. For example, in the case of equal harmonic traps of 
strength Q, such that = fi (i.e., for soliton motion near the trap center), the 

oscillation frequency of the dark-bright soliton resulting from Eq. (11321) reads, 

\ 1/2 

Q osc = ^=\l ; 2±= = xo I • (133) 




4v//i+(iV b /4) 2 

It is clear that Eq. fl 1 3 3 [) shows that the oscillation frequency is down-shifted as compared 
to the characteristic value of f2/\/2, i-e., the dark-bright pair executes slower oscillations, 
as the bright component is enhanced. 
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The predictions of Ref. [131] can directly be compared to the findings of a Hamburg 
experiment [67J, where long-lived dark-bright matter- wave solitons were observed in a 
two-component quasi-lD 87 Rb BEC. In particular, using the phase-imprinting method, 
a dark soliton was created in one spin state of the BEC and the density dip was filled 
by atoms, forming the bright soliton, in another spin state of the BEC (note that the 
number of atoms of the bright soliton was ~ 10% of the total number of atoms). 
The created dark-bright soliton was then observed to perform slow oscillations with a 
frequency 0.24S7, which is quite smaller than the frequency of the corresponding single 
dark soliton in the same setting. Moreover, due to the initial state preparation, an 
extra dark soliton was generated, which was allowed to interact with the co-existing 
dark-bright soliton; it was observed that this individual dark soliton was reflected off 
the slower dark-bright one, with the process resembling a hard-wall reflection. 

6.1.2. Dark solitons in spinor condensates. The spin degree of freedom of spinor BECs 
gives rise to important new phenomena (including, among others, the formation of spin 
domains [242 j . spin textures [282] and vortices [283] . as well as spin oscillations [284] ). 
which are not present in other types of BECs. Generally, a spinor condensate formed 
by atoms with spin F can be described in the framework of mean-field theory by a 
(2F + l)-component macroscopic wave function; accordingly, a spinor F = 1 condensate 
is characterized by a vector order parameter, with the three components corresponding 
to the values of the vertical spin projection, = —1,0,-1-1. In a quasi-lD setting, 
the pertinent system of coupled GP equations for the wave functions ip±ifi(z,t) can be 
expressed in the following dimensionless form (see, e.g., Refs. |123[l285jl286] ). 

i$V±i = «V±i + W±i\ 2 + IV'oI 2 - |<M 2 )^±i + #0^1, (134) 
id t ^o =^o + W-i| 2 + |V>+i| 2 )V>o + (135) 
Here, H. = — (l/2)d 2 + (l/2)Q 2 z 2 + n, with Q being the normalized trap strength 
and n = \ip-i\ 2 + | V^o 1 2 + iV'+il 2 the total density, while 5 = (a 2 — a )/(a + 2a 2 ) 
where ao and a 2 are the s-wave scattering lengths in the symmetric channels with total 
spin of the colliding atoms F = and F = 2, respectively. Actually, the parameter 
5 represents the ratio of the strengths of the spin-dependent and spin-independent 
interatomic interactions, and may take negative or positive values for ferromagnetic 
or anti-ferromagnetic (alias polar) spinor BECs, respectively. Typically, in the relevant 
cases of 87 Rb and 23 Na atoms with F = 1, 5 = -4.66 x 10~ 3 [2E7] and 5 = +3.14 x 10~ 2 
[288J; nevertheless, the above values may in principle be modified by employing the 
so-called confinement- induced Feshbach resonance [289] . 

In the limiting case of 5 = (and in the absence of the potential), the system of 
Eqs. (]134p — (I135p is reduced to the completely integrable Manakov system. On the other 
hand, as shown in Ref. [290] , another completely integrable version of Eqs. f ll34p -( fl~35i) 
corresponds to the case 5 = 1 (i.e., for interatomic and anti- ferromagnetic interactions of 
equal magnitude): in this case, the resulting matrix NLS equation with non- vanishing 
boundary conditions is completely integrable by means of the 1ST method [291] and 
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admits exact analytical vector iV-dark soliton solutions (i.e., single- and multiple- vector 
dark solitons of the dark-dark-dark type in terms of the mp = —1,0,-1-1 spinor 
components) |292] (see also Ref. p93j). The one-dark soliton state of this system can 
be classified as (a) ferromagnetic (i.e., with nonzero total spin), which has domain- 
wall shaped wave functions, and (b) polar (i.e., with zero total spin), characterized by 
the familiar hole soliton profile. Note that the collisions of two solitons give rise to 
interesting spin-dependent phenomena, such as spin- mixing or spin-transfer [292J. 

In the physically relevant case of small 5, mixed dark-bright solitons of the dark- 
dark-bright or bright-bright-dark type (again in terms of the nip = —1,0,-1-1 spinor 
components) were also predicted to occur in anti-feromagnetic spinor F — 1 BECs [123j . 
In the small-ampliude limit (and in the absence of the trap), these solitons were found 
to obey the completely integrable Yajima-Oikawa system p94j, by means of which it 
was found that the functional form of the dark and bright components is similar to the 
one in Eqs. fll30p - fll31l) . Numerical simulations in Ref. |123] demonstrated that, for 
small-amplitudes, such dark-bright solitons feature genuine soliton behavior (i.e., they 
propagate undistorted and undergo quasi-elastic collisions), while for moderate and large 
amplitudes (and also for large values of 8) they can exist as long-lived objects as well. 
Furthermore, for sufficiently small number of atoms of the bright soliton, the bright 
component (s) are guided by the dark one(s), and the vector soliton performs harmonic 
oscillations; the oscillation frequency is different for small- and moderate-amplitude 
solitons, and it is respectively given by: 

Wosc = ^=(1 _ "ov^) - e , tu osc = fi osc (l - «i5 2 ) + e x . (136) 

In the above expression, Q osc is given by Eq. (11331) . while the constants ao,i and e 0; i 
(with eo,i <C «o,i) depend on the normalized number of atoms of the bright component. 
It is clear that the characteristic oscillation frequency of the dark soliton (VL/y/2), as 
well as the result of Ref. [131] . is modified by the spin- dependent interactions that are 
present in the case of a spinor F = 1 BEC. 

6.2. Matter-wave interference and dark solitons. 

Matter- wave interference experiments (see, e.g., the seminal work of Ref. |295j ) are 
known to demonstrate, apart from self-interference, the interference between two BECs 
confined in a trap, divided into separate parts by means of a barrier potential induced by 
a laser beam. In particular, the BECs are left to expand and overlap forming interference 
fringes, similar to the ones known in optics. Much interest has been drawn to a better 
understanding of this fundamental phenomenon, especially as concerns the coherence 
properties of the interfering BECs. In that regard, it is worth mentioning that the 
(incorrect) assumption that "when the interfering BECs have fixed atom numbers, 
there can be no phase", was resolved - - shortly after the experimental realization 
of BECs [361438] - in Ref. [296] . On the other hand, since most of the relevant 
experimental findings can be quantitatively reproduced in the framework of the GP 
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mean- field theory [297J, we will proceed by adopting this approach in order to discuss 
the connection between dark solitons and matter- wave interference. 

An interesting variation of the interference process, which is naturally attributed 
to the inherent nonlinearity of BECs due to interatomic interactions, is that - 
under certain conditions - - the collision of two initially separated condensates can 
lead to the creation of dark solitons. This "nonlinear interference" effect was first 
observed in simulations |139] , and was subsequently analyzed theoretically [298J . Other 
studies, basically relying on the self-interference of BECs, have also been proposed 
as well [220T [223jl299j . Importantly, relevant recent experiments employing this, so- 
called, matter-wave interference method have already been reported, demonstrating the 
generation of vortices [300J (see also theoretical work in Refs. |223[l301|l302j ) and dark 
solitons [691472] (see also the experiment of Ref. [62]). 

To get a deeper insight into the physics of the matter-wave interference process, let 
us follow the arguments of Ref. [298] and consider the interference between two separated 
quasi- ID BECs colliding in the presence of a harmonic trap. There exist two different 
regimes characterizing this process, namely a linear and a nonlinear one, depending on 
the competition between the kinetic and the interaction energies. In the linear regime, 
the kinetic energy of the condensates exceeds the nonlinear interaction energy of the 
atoms. In this case, and at any time t, the total wave function of the system can be 
well approximated by a linear superposition of the wave functions that each individual 
condensate would have at t. The two initially well-separated BECs interfere at the trap 
center, produce a linear interference pattern, and then separate again regaining their 
initial shape. The fringe spacing I of the interference pattern is determined by the 
k vector that each individual condensate would have if performing a dipole oscillation 
alone in the trap and, at the time of maximum overlap, I = it/ D(h/2mu z ) (here D is the 
initial distance between the condensates). It is clear that the higher the kinetic energy 
is, the higher the number of fringes and the smaller the fringe spacing is. Approximating 
the individual wave functions in the TF limit, it can be found that the kinetic energy 
(estimated from the curvature of a cos 2 interference pattern) exceeds the peak nonlinear 
energy (at the center of the fringes) when the initial distance between the two BECs 
exceeds critical distance, namely D > D c , where D c is given by [298J: 



Here, N is the number of atoms, a the s-wave scattering length, oo z the longitudinal 
trap frequency, and m the atomic mass. If the above condition is not fulfilled, namely 
D < D c , then the system enters in the nonlinear regime. In the latter, the interference 
pattern consists of stable fringes with a phase jump of order of 7r across them, which can 
naturally be identified as genuine dark solitons. Notice that in the nonlinear regime, 
the initially individual condensates instead of reforming as separate objects, they form 
a combined condensate undergoing a quadrupole oscillation. 

The results of Heidelberg experiments [69|I7T] can be compared directly to the above 
theoretical predictions. In these experiments, dark solitons were created by releasing a 
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87 Rb BEC from a double-well trap into a harmonic trap in the dimensionality crossover 
regime from ID to 3D. For the parameters used, the initial distance of the individual 
condensates was approximately five times smaller than the critical distance and the trap 
frequencies were ramped, with ramping times chosen so as to minimize the excitation of 
the quadrupole mode. It was shown that, in accordance to the observations of Ref. [62J, 
the number of created solitons is even for a zero phase- difference between the two initially 
separated condensates, while it is odd for a phase-difference close to n. If the phase- 
difference is exactly equal to it, a standing (black) dark soliton in the middle of the trap 
is always created. Notice that the total number of the created solitons depends on the 
momentum of the merging condensates, which may be controlled by varying, e.g., the 
distance between the condensate fragments, the number of atoms, or the aspect ratio of 
the trap (SUEDES]. 

6.3. BEC superfluidity and dark solitons 

A flow past an obstacle is known to be one of the most fundamental contexts for studying 
superfluidity. Particularly, according to the Landau criterion for superfluidity |303j . a 
superfluid flow past an obstacle, is stable (unstable) for group velocities smaller (larger) 
than the speed of sound. Actually, breakdown of superfluidity is caused by the opening 
of channels for emission of excitations in the fluid, whose formation manifests itself 
as an effective dissipation. In the BEC context, early experiments from the MIT 
group |187[l304j demonstrated the onset of dissipation induced by the motion of an 
obstacle (in the form of a strongly repulsive dipole beam). From a theoretical standpoint, 
the problem can be studied by using a NLS (or a GP) equation that includes a localized 
external potential of the form V(r — vt) (with V(r) — » as |r| — » oo) accounting for 
the presence of the obstacle moving with velocity v; this potential may be naturally 
superimposed to the usual trapping potential confining the condensate. In relevant 
earlier studies [305J, where the NLS equation as a model of superflow was used, vortex 
formation induced by the superfluid flow around an obstacle was predicted. 

The lower-dimensional setting, namely the ID flow of a repulsive NLS fluid in the 
presence of an obstacle, was also studied (306[l307j (see also Ref. |308j ). Specifically, in 
Ref. [306] it was shown that below an obstacle-dependent critical velocity, there exists a 
steady dissipationless flow solution, which disappears at the critical velocity by merging 
with an unstable solution in a saddle-node bifurcation. This unstable solution represents 
the transition state for emission of dark solitons, which are repeatedly generated above 
the critical velocity. In fact, the onset of dissipation corresponds to nonstationary flow 
with a wake asymptotically extending upstream to infinity, and downstream periodic 
emission of dark solitons [307J. Note that in both Refs. [306J and [307] the critical 
velocity was found to be smaller than the speed of sound, a result that may be explained 
by the fact that, in the region of the potential, the local fluid velocity can reach values 
higher than the local sound velocity (critical velocity values smaller than the speed of 
sound were also observed in Refs. |187[I304] ). 
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The above results paved the way for a better understanding of the BEC flow past 
an obstacle and inspired further investigations [3091- 131 1] . Importantly, in an recent 
experiment |66j, the BEC flow induced by a broad, penetrable barrier (in the form of 
a laser beam) swept through an elongated 87 Rb condensate was systematically studied: 
it was demonstrated that at slow barrier speeds the flow is stable, at intermediate 
speeds becomes unstable and dark soliton generation is observed, while at faster speeds, 
remarkably, soliton formation completely ceases. Both repulsive and attractive barriers 
were used in the experiment of Ref . [66] , and were found to lead to dark soliton formation; 
additionally, it was also found that the critical velocity for the breakdown of the BEC 
superfluidity and soliton generation was smaller than the speed of sound. Note that in a 
recent work [312] . velocity regimes similar to the ones found in Ref. (66] were analytically 
predicted by using a hydrodynamic approach. 

As shown theoretically [306, 307] I309H311] and demonstrated experimentally [66] . 
dark solitons (and vortices) are formed if the size of the "hypersonic" obstacle is of 
the order of, or greater than, the characteristic healing length of the condensate. On 
the other hand, if the size of the obstacle is much smaller than the healing length, the 
main loss channel, which opens at supersonic velocities of the obstacle, corresponds 
to the Cerenkov emission of Bogoliubov's excitations |313] . Notice that in the case of 
large hypersonic obstacles, two dispersive shock waves, which start propagating from 
the front and the rear parts of the obstacle, are formed. Far from the obstacle, the 
shock front gradually transforms into a linear "ship wave" located outside the Mach 
cone [314H317] . whereas the rear zone of the shock is converted into a "fan" of oblique 
dark solitons located inside the Mach cone 13 1811320] (see also relevant experimental 
results in Refs. |314U321j ). An important result reported in Ref. |322j is that although 
such dark solitons are unstable in higher-dimensional settings with respect to transverse 
perturbations (see Sec. I5.ip . the instability becomes convective — rather than being 
absolute — for sufficiently large flow velocities and, thus, dark solitons are effectively 
stable in the region around the obstacle. 

The flow of a multi-component BEC past an obstacle was also studied, and the cases 
of a two-component [270|,l271[ l323j and a spinor F — 1 condensate [324J were analyzed. 
It is interesting to note that, as shown in Ref. |270] in the case of a two-component BEC, 
the existence of two different speeds of sound provides the possibility for three dynamical 
regimes: when both components are subcritical, nucleation of coherent structures does 
not occur; when both components are supercritical they both form dark solitons in ID 
and vortices or rotating vortex dipoles in 2D; in the intermediate regime, the nucleation 
of a dark-anti-dark soliton in ID or a vortex-lump configuration in 2D is observed. 
Furthermore, as shown in Ref. [271] , dark solitons can be convectively stabilized in the 
2D setting at sufficiently high values of the obstacle velocity, similarly to the case of 
one-component BECs [322J. 
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6-4- Matter-wave dark solitons in optical lattices. 

Bose-Einstein condensates loaded into periodic optical potentials, so-called optical 
lattices (OLs), have attracted much attention as they demonstrate rich physical 
properties and nonlinear dynamics (see, e.g., Refs. [43,83,325f l328j for reviews). Optical 
lattices are generated by a pair of laser beams forming a standing wave which induces a 
periodic potential; thus, for a BEC confined in an optical lattice, the trapping potential 
in the GP model can be regarded as a superposition of a harmonic (magnetic or optical) 
trap and a periodic potential. Particularly, in a quasi- ID setting (generalization to 
higher- dimensional cases is straightforward) - - cf. Eq. (17H|) - - the trap takes the 
following dimensionless form (see, e.g., Ref. |44] ) : 

V(z) = ~OV + V Q cos 2 {kz), (138) 

Here, Q and V denote, respectively, the harmonic trap and OL strengths, L = n/k = 
(A/2) sin (0/2) is the periodicity of the lattice, with A being the common wavelength of 
the two interfering laser beams, and the angle between them. In some cases (as, e.g., 
in the experiments of Refs. [32911330] ). the harmonic potential is very weak as compared 
to the optical lattice and, thus, it can be ignored. Then, the stationary states of the 
pertinent GP equation — including solely the OL potential — can be found in the 
form of infinitely extended waves, with the periodicity of the OL, known as nonlinear 
Bloch waves (see, e.g., Ch. 6 in Ref. [13] and references therein). In the same case 
(i.e., in the absence of the harmonic potential), if the OL is very deep (compared to 
the chemical potential), the strongly spatially localized wave functions at the lattice 
sites are approximated by Wannier functions (see, e.g., Ref. |331j ) and the tight-binding 
approximation can be applied; then, the continuous GP equation is reduced to the 
discrete NLS (DNLS) equation (see, e.g., Refs. (BHS31E25] and Ch. 13 in Ref. pj, as 
well as Refs. |332[ [333j for reviews for the DNLS model). Dark solitons, which may 
naturally exist in all of the above settings and combinations thereof, have been studied 
both in combined harmonic and OL potentials |183|.ll84|,l334|,l335j and in optical lattices 
(in the absence of the harmonic trap). In the latter case, various studies have been 
performed in the frameworks of the continuous GP equation, as well as its tight-binding 
approximation counterpart [336H340] . Notice that matter- wave dark solitons have also 
been studied in double-periodic optical superlattices f [340U341] . while there exists a vast 
amount of work concerning dark solitons in periodic media arising in various contexts, 
such as nonlinear optics [343-348J, solid-state physics |349j and the theory of nonlinear 
waves [350J. 

6.4-1- Dark solitons in combined harmonic and OL potentials. The stability of matter- 
wave dark solitons in the combined harmonic and OL potential, was first studied in 
Ref. [334] by a means of a BdG analysis that was performed in the framework of both 

f Such a potential has the form V(z) = Vi cos{k\z) + V2 cos^z); where ki and k% > fci are the primary 
and secondary lattice wavenumbers, and V\ and V2 are the associated sublattice amplitudes [342) . 
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the continuous quasi-lD GP equation and its DNLS counterpart. It was found that in 
the discrete model stationary dark solitons located at the minimum of the harmonic 
trap are, generally, subject to a weak oscillatory instability, which manifests itself as a 
shift of the soliton from its initial location, accompanied by quasi-periodic oscillations. 
On the other hand, in the continuous GP model, dark solitons may be stable, with the 
(in) stability determined by the period and amplitude of the OL. In any case, the dark 
solitons are robust and if the oscillatory instability is present, it sets in at large times. 

The dynamics of dark solitons in the combined harmonic and OL potential can be 
studied upon distinguishing physically relevant cases, depending on the competition of 
the characteristic spatial scales of the problem [184J. Particularly, assuming that the 
harmonic trap varies slowly on the soliton scale, i.e., w — 1/ cos</> ~ £ <C Q^ 1 (where 
w is the soliton width for chemical potential /i = 1, <fi is the soliton phase angle, and 
£ the healing length), the following three cases can readily be identified: (a) the case 
of a long-period OL, with L ^> £, (b) the case of a short-period OL, with L <C £, 
and (c) the intermediate case, with L ~ £. Then, if the OL strength is sufficiently 
small, the soliton dynamics in cases (a) and (b) can be treated in the framework of the 
adiabatic approximation (see Sec 14. 2ft . Particularly, as shown in Ref. |184j . case (a) can 
be studied by means of the Hamiltonian approach of the perturbation theory, and case 
(b) by means of a multi-scale expansion method (treating k~ x as a small parameter); this 
way, it can be shown that in both cases the dark soliton behaves as an effective classical 
particle, performing harmonic oscillations in the presence of the trap of Eq. (11381) . The 
oscillation frequency, which is different from its characteristic value Q/y2, is modified 
by the presence of the lattice according to the equations: 



for cases (a) and (b), respectively. As concerns the more interesting case (c) (see 
Refs. |183[I184] ). it can be shown that if the dark soliton is initially placed quite close 
to the bottom of a well of the OL potential, it remains there for a rather long time; 
eventually, however, it escapes due to the radiation-loss mechanism, and then performs 
large-amplitude oscillations in the condensate. Furthermore, if the harmonic trap is 
weak enough, the soliton eventually decays. In fact, as discussed in Ref. |183j . the OL 
causes a dynamical instability (because the dark soliton has to "traverse" the potential 
humps caused by the lattice) resulting in a faster decay of the soliton than if it was 
evolving in the presence of the harmonic trap only: the presence of the lattice dephases 
the sound waves emitted by the soliton, hence reducing the effectiveness of the soliton to 
get stabilized by reabsorbing the sound waves (see discussion in Sec. I4.4p . Nevertheless, 
according to the observations of Ref. [184] that the soliton may remain stationary for 
a relatively long time, in Ref. |335] (see also Ref. |184j ) it was proposed that a time- 
dependent OLs may either (i) capture a moving dark soliton or (ii) capture and drag a 
stationary soliton, bringing it to a pre-selected final destination. Notice that the transfer 
mechanism is robust as long as adiabaticity of the process is ensured (i.e., for sufficiently 
small speeds of the moving OL). 
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6.4-2. Dark solitons in optical lattices and superlattices. As mentioned above, dark 
solitons in OLs and superlattices have also been studied, in the absence of the harmonic 
trap, in the frameworks of the continuous and discrete NLS equations. Particularly, in 
Ref. [336], a DNLS model was derived in the tight-binding approximation — i.e., for a 
single-isolated band in the Floquet-Bloch spectrum — which was used to study matter- 
wave dark solitons. Later, in Ref. [337] , a continuous coupled-mode model was used to 
study the existence and stability of, so-called, dark lattice solitons, while a more general 
analysis was presented in Ref. [338] : in that work, a continuous GP model with periodic 
potential was shown to support stable stationary dark solitons for both attractive and 
repulsive interatomic interactions, which were found numerically |338] (see also relevant 
results in |339j ). 

In a more recent work |340j . where both regular optical lattices and superlattices 
were considered, it was shown that each type of nonlinear Bloch wave can serve as 
a stable background supporting dark solitons. This way, different families of dark 
solitons, originating within the bands of the Floquet-Bloch spectrum, were found and 
their dynamical properties were analyzed. In particular, considering the continuous 
analogue of the Peierls-Nabarro potential (see, e.g., Ref. |351j ) in discrete lattices, it was 
shown that the mobility and interaction properties of the dark solitons can be effectively 
controlled by changing the structure of the optical superlattice; moreover, following the 
ideas of Ref. [335], time-dependent superlattices were also shown to control the static 
and dynamical properties of matter- wave solitons [341J. 

Here we should point out that all the above mentioned studies on the dynamics 
of matter-wave dark solitons in optical lattices were carried out in the framework 
of the GP mean-field theory. Nevertheless, it is worth emphasizing that the GP 
equation is inadequate for dealing with several important aspects of ultracold bosons 
in optical lattices, such as the superfluid-to-Mott insulator phase transition (see, e.g., 
Refs. [352(1353] ) or, more generally, strong correlation effects (see, e.g., the review [328] ). 
Thus, more recently, studies on the quantum dynamics of dark solitons have started 
to appear. In that regard, it is relevant to mention that matter-wave dark solitons 
were studied in the context of the Bose- Hubbard model [354[ l355j. and it was found 
that dark soliton collisions become inelastic, in strong contrast to the predictions of 
mean-field theory. A conclusion of the above works [354U355] is that the lifetime and 
collision properties of matter-wave dark solitons in optical lattices may provide a clear 
signature of quantum effects. Additionally, in another recent work [356] . dark solitons 
were studied in the vicinity of the superfluid-to-Mott insulator transition; particularly, 
in this work [356] . antisymmetric eigenstates corresponding to standing solitons, as well 
as propagating solitons created by phase-imprinting, were presented and the soliton 
characteristics were found to depend on quantum fluctuations. 

From the viewpoint of experiments, trains of stationary dark solitons were observed 
in a 87 Rb condensate confined in a 3D harmonic trap and a ID OL [M]. The underlying 
mechanism for the formation of such structures were multiple Bragg reflections caused 
by displacing the harmonic trap and, thus, setting the BEC into motion. Due to the 
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dimensionality of the system, the solitons were found to be subject to the snaking 
instability, giving rise to the subsequent formation of vortex rings (see Sec. I5.1[) . similarly 
to the observations of the pertinent JILA experiment (but without the OL) [48j. 



6. 5. Matter-wave dark solitons at finite temperatures. 

So far, we have considered the stability and dynamics of matter-wave dark solitons at 
zero temperature, T = 0. Nevertheless, as experiments are obviously performed at finite 
temperatures, it is relevant to consider the dissipative instability of dark solitons induced 
by the thermal excitations that naturally occur. This problem was first addressed in 
Ref. [112] . where a kinetic-equation approach, together with a study of the Bogoliubov- 
de Gennes (BdG) equations, was used. In this work, it was found that the dark soliton 
center obeys an equation of motion which includes an anti- damping term [similarly to 
Eq. fll05p ], which is nonzero (zero) for finite (zero)-temperature. The behavior of the 
solutions of this equation of motion incorporating the anti-damping term can be used 
to explain — at least qualitatively — the soliton dynamics observed in experiments: 
solitons either decay fast (for high temperatures) [I5],ll6lll9] or perform oscillations (for 
low temperatures) [6TW69| l7T] of growing amplitude and eventually decay, so that the 
system finally relaxes to its ground state (see also discussion below for the role of the 
anti-damping term). 

Dark soliton dynamics in BECs at finite temperatures was also studied in other 
works by means of different approaches. In particular, in Ref. [185] . the problem 
was treated in the framework of a mean-field model, namely the so-called dissipative 
GP equation (see below); this equation incorporates a damping term, first introduced 
phenomenologically [357] and later justified from a microscopic perspective (see, e.g., 
the review [358] ). On the other hand, in Refs. [359(1360] the same problem was studied 
numerically, using coupled Gross-Pitaevskii and quantum Boltzmann equations, which 
include the mean field coupling and particle exchange between the condensate and 
the thermal cloud. Furthermore, in Refs. [5811361] . finite-temperature dynamics of 
dark solitons was studied by means of the so-called stochastic GP equation (see, e.g., 
Ref. [358] ). while in Ref. [362] quantum effects on dark solitons were additionally studied 
in the framework of the truncated Wigner approximation (see, e.g., Refs. [363,364| for 
this approach); we also note that in the recent work [365J, the dissipative dynamics of a 
dark soliton at temperatures T, lower than the chemical potential /i of the background 
Bose liquid, was studied. In the work [361J, it was shown that for sufficiently low 
temperatures and certain parameter regimes, averaged dark soliton trajectories obtained 
by the stochastic GP equation, are in a very good agreement with results obtained by 
the dissipative GP model. Thus, the results of Ref. [361] indicate that the use of the 
dissipative GP equation in studies of dark solitons in finite-temperature BECs (a) can 
reasonably be justified from a microscopic perspective, and (b) allows for an analytical 
description of the problem, by employing techniques exposed in Sec. 14.21 provided that 
the dynamics of the thermal cloud does not play the dominant role. 
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To be more specific, we follow Ref. [361] . and express, at first, the dissipative GP 
model in the following dimensionless form, 

'1 



f 2 z + V(z) + \^-fi 



1>, 



(140) 



where units are the same to the ones used for Eq. f[24|) . and the dimensionless 
dissipation parameter 7 can be connected with temperature by means of the relation 
7 oc (ma 2 kBT)/(TTh 2 ), where ks is Boltzmann's constant. Dark matter-wave soliton 
dynamics can be studied analytically in the framework of Eq. (I140p . by employing the 
Hamiltonian approach of the perturbation theory for dark solitons (see Sec. 14.2.21) . In 
particular, we assume that the condensate dynamics involves a fast scale of relaxation of 
the background to the ground state, while the dark soliton subsequently evolves on top 
of the relaxed ground state; then, it is possible to derive the perturbed NLS Eq. (1541) for 
the dark soliton wave function ip s , with a perturbation Q{jp s ) [cf. Eq. (155]) ] incorporating 
the additional term 2 r yfid t ip s . Then, following the procedure of Sec. 14.2.21 we end up 
with the following equation of motion for the dark soliton center z 
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In the case of nearly-black solitons with dz$/dt sufficiently small, Eq. (I14ip can be 
reduced to the following linearized form (similar to the equation of motion of Ref. |112j ): 

z = 0. (142) 



d 2 z 



dt 2 



2 dz , 



7^) 



In the limiting case of zero temperature, 7 = 0, Eq. (I142j) is reduced to Eq. (1901) [for 
the harmonic trap V(z) = (l/2)Q 2 z 2 ]. On the other hand, at finite temperatures, 
7^0, Eq. (I142p incorporates the anti-damping term (oc —dzo/dt). Although it may 
sound counter-intuitive, this term describes the dissipation of the dark soliton due to 
the interaction with the thermal cloud: in fact, this term results in the acceleration of 
the soliton towards the velocity of sound, i.e., the soliton becomes continuously grayer 
and, eventually, the soliton state transforms to the ground state of the condensate. 

Explicit solutions of Eq. ( I142p can readily be obtained in the form of z oc exp(s 12 t), 
where Si j2 are the roots of the auxiliary equation s 2 — (2/3)jijls + (Q/\/2) 2 = and are 
given by: 
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In Ref. |361j (see also relevant work in Ref. |366] ) , the temperature dependence of these 
eigenvalues associated with the dark soliton dynamics was compared to the temperature 
dependence of the eigenvalues of the pertinent anomalous mode of the system |361j . and 
the agreement between the two was found to be excellent. Both the motion eigenvalues 
[cf. Eq. ( I143p ] and the anomalous mode eigenvalues (derived by a BdG analysis) 
undergo Hopf bifurcations as the dissipation (temperature) is increased/decreased - 
leading to an exponential/oscillatory instability of the dark soliton — with the respective 
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bifurcation diagrams being almost identical. A similar situation occurs in the case of 
multiple-dark solitons as well: as shown in Ref. [366], eigenvalues derived by coupled 
equations of motions [similar to the ones in Eqs. (1122^ — (I123p ] for two- or three-dark 
solitons, were again found to be almost identical to the ones of the anomalous modes of 
the system. 



7. Conclusions and perspectives 

We have presented the recent progress on the study of dark solitons in atomic BECs, 
including analytical, numerical and experimental results. In fact, although the main 
body of this work was basically devoted to the theoretical aspects of this topic, we have 
tried to connect the theoretical results to pertinent experimental observations. In that 
regard, we have particularly tried to highlight the close connection between theory and 
experiments and the reasonable agreement between the two. 

Matter- wave dark solitons were predicted to occur in BECs as early as 1971 [9], 
but were observed in experiments only 28 years later, in 1999 [15]. Although, till then, 
dark solitons had already a relatively long history in the context of nonlinear optics 
(where they were first observed in experiments on 1987 pQ and studied extensively 
in theory during the next years [31]), one can readily realize an emerging interest in 
them: during the last decade, there have been more than ten experiments on dark 
solitons [^5l - H9ll63ti7T] . and half of them have been conducted very recently [661471] . with 
an unprecedented control over both the condensate and the solitons. As the experimental 
developments continuously inspire — and, at the same time, are guided by — a huge 
number of relevant theoretical works, one may expect that the interest in matter-wave 
dark solitons will still be growing in the near future. 

Although there has been a tremendous progress on our understanding of matter- 
wave dark solitons in atomic BECs over the last years, many important issues remain 
to be addressed or studied in a more systematic way. A relevant list is appended below. 

• Beyond Mean-Field. Matter-wave dark solitons, being fundamental nonlinear 
macroscopic excitations of BECs, play an important role at probing the properties of 
the condensates at the mesoscale (see, e.g., discussion in Ref. |59j). In that regard, 
a quite interesting research direction is the study of these nonlinear structures, 
both in theory and in experiments, in various settings and regimes where thermal 
and quantum effects are important. In fact, mean- field theory can only account 
for averaged results (e.g., soliton decay times [36 lj ). whereas recent experiments 
[6711691171] indicate shot-to-shot variations that could be accounted for by stochastic 
approaches. There exist various experimentally relevant settings — such as the 
ones where the number of atoms is small, particularly those in optical lattices or at 
very low temperatures — which enhance the importance of quantum fluctuations; 
therefore, the latter should be appropriately included. One interesting question 
concerns, for example, the issue of the filling of the dark soliton due to averaging 
based on thermal or quantum fluctuations [169113671 1368J and its relation to the 
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measurement process [369] . Extending this argument, one could use dark soliton 
experiments to test the regimes of validity of conventional mean-field theories, a very 
interesting and fundamental topic in its own right. From a theoretical standpoint, 
the above directions seem to be a natural next step in the study of BECs and their 
excitations; in fact, relevant work — based on various approaches beyond the mean- 
field approximation — has already started (see, e.g., Refs. |354H356[I361[I362[I365] ) 
and is expected to continue even more intensively in the near future. 

• Mathematical analysis. Even in the framework of the mean-field approximation, 
there exist several theoretical problems which remain unsolved or should be 
investigated in more detail. A pertinent example is the study of the persistence 
and stability of dark solitons in the presence of confining or periodic potentials: 
as mentioned in Sec. 14.51 rigorous results have only been obtained for small, 
bounded and decaying potentials [161], while an analysis of other cases is still 
missing. Furthermore, there is still work to be done as concerns the development 
of perturbation theories for multiple dark solitons, for dark solitons in multi- 
component systems, dissipative systems, and others. 

• Further experiments. From the viewpoint of experiments, the recent observations 
of long-lived matter- wave dark solitons [67H691 171] suggest many other possible 
experimental investigations. In fact, there are many interesting problems related to 
dark solitons, which require experimental studies. These include (a) the influence 
of thermal and quantum fluctuations on dark solitons (as indicated above), (b) 
investigation of states composed by a large number of dark solitons (including, 
so-called, "soliton gases" - - see, e.g., Refs. |370[l371j ). (c) observation of vector 
solitons, such as dark-dark and dark-ant i- dark solitons in two-component BECs 
or vector solitons with at least one component being a dark soliton in spinor 
BECs [T23U265U267U292U293] . (d) interactions of dark solitons with potential barriers 
and studies of the reflectivity /transmittivity of dark solitons [731 11371 [1461 1159] . 
(e) manipulation of dark solitons in collisionally inhomogeneous environments 
|188U189j . or by means of time-dependent optical lattices [184l [335 ,341]. and others. 

• Applications. Apart from basic theory and relevant experiments, an important 
question concerns possible applications of matter-wave dark solitons. Although 
there exist some works indicating the importance of dark solitons in atomic matter- 
wave interferometers in the nonlinear regime [60-63, 299J - a direction which is 
expected to further be explored - - other potential applications (similar to the 
ones related to optical dark solitons [SI EH]) remain to be investigated. As an 
example we note that matter-wave bright-dark vector solitons (which have already 
been observed [67]) in pseudo-spinor or spinor BECs may provide the possibility 
of all-matter-wave waveguiding: in such a situation, the dark soliton component 
could build an effective conduit for the bright component, similar to the all-optical 
waveguiding proposed in nonlinear optics [31]. Waveguides of this kind would be 
useful for applications, such as quantum switches and splitters emulating their 
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optical counterparts |372j . 

• Ultracold Fermi gases. We finally note that, so far, matter-wave dark solitons 
have mainly been studied in the context of ultra-cold Bose gases. Nevertheless, 
recent progress in the area of ultra-cold Fermi gases (see, e.g., Refs. [37311374] for 
recent reviews), suggest that (similarly to vortices) dark solitons may be relevant 
in this context as well. In fact, pertinent theoretical studies have already started 
to appear |375H377j . but there is still much work to be done towards this direction, 
both in theory and in experiments. 
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